Bayesian GLMM Part2

Author

Murray Logan

Published

07/07/2025

1 Preparations

Load the necessary libraries

library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(standist) # for exploring distributions
library(HDInterval) # for HPD intervals
library(posterior) # for posterior draws
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(DHARMa) # for residual diagnostics
library(tidybayes) # for more tidying outputs
library(ggeffects) # for partial plots
library(broom.mixed) # for tidying MCMC outputs
library(patchwork) # for multiple plots
library(ggridges) # for ridge plots
library(bayestestR) # for ROPE
library(see) # for some plots
library(easystats) # framework for stats, modelling and visualisation
library(modelsummary)
source("helperFunctions.R")

2 Scenario

To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.

A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.

**Figure** 2: Sampling design for the norin data set

Table 1: Format of norin.csv data files

FISHID	MASS	TRIAL	SMR_contr	CHANGE
1	35.69	LowSalinity	5.85	-31.92
2	33.84	LowSalinity	6.53	2.52
3	37.78	LowSalinity	5.66	-6.28
..	..	..	..	..
1	36.80	HighTemperature	5.85	18.32
2	34.98	HighTemperature	6.53	19.06
3	38.38	HighTemperature	5.66	19.03
..	..	..	..	..
1	45.06	Hypoxia	5.85	-18.61
2	43.51	Hypoxia	6.53	-5.37
3	45.11	Hypoxia	5.66	-13.95

Table 2: Description of the variables in the norin data file

FISHID	Categorical listing of the individual fish that are repeatedly sampled
MASS	Mass (g) of barramundi. Covariate in analysis
TRIAL	Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen.
SMR_contr	Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia)
CHANGE	Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'.

3 Read in the data

norin <- read_csv("../data/norin.csv", trim_ws = TRUE)

Rows: 180 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): TRIAL
dbl (4): FISHID, MASS, SMR_contr, CHANGE

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(norin)

Rows: 180
Columns: 5
$ FISHID    <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ MASS      <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.9…
$ TRIAL     <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity", …
$ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589, …
$ CHANGE    <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.0…

## Explore the first 6 rows of the data
head(norin)

# A tibble: 6 × 5
  FISHID  MASS TRIAL       SMR_contr CHANGE
   <dbl> <dbl> <chr>           <dbl>  <dbl>
1      1  35.7 LowSalinity      5.85 -31.9 
2      2  33.8 LowSalinity      6.53   2.52
3      3  37.8 LowSalinity      5.66  -6.28
4      4  26.6 LowSalinity      6.28  -4.35
5      5  37.6 LowSalinity      4.41  -3.07
6      6  37.7 LowSalinity      4.82 -15.1

str(norin)

spc_tbl_ [180 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ FISHID   : num [1:180] 1 2 3 4 5 6 7 8 9 10 ...
 $ MASS     : num [1:180] 35.7 33.8 37.8 26.6 37.6 ...
 $ TRIAL    : chr [1:180] "LowSalinity" "LowSalinity" "LowSalinity" "LowSalinity" ...
 $ SMR_contr: num [1:180] 5.85 6.53 5.66 6.28 4.41 ...
 $ CHANGE   : num [1:180] -31.92 2.52 -6.28 -4.35 -3.07 ...
 - attr(*, "spec")=
  .. cols(
  ..   FISHID = col_double(),
  ..   MASS = col_double(),
  ..   TRIAL = col_character(),
  ..   SMR_contr = col_double(),
  ..   CHANGE = col_double()
  .. )
 - attr(*, "problems")=<externalptr>

norin |> datawizard::data_codebook()

norin (180 rows and 5 variables, 5 shown)

ID | Name      | Type      | Missings |           Values |          N
---+-----------+-----------+----------+------------------+-----------
1  | FISHID    | numeric   | 0 (0.0%) |          [1, 60] |        180
---+-----------+-----------+----------+------------------+-----------
2  | MASS      | numeric   | 0 (0.0%) |   [22.52, 64.07] |        180
---+-----------+-----------+----------+------------------+-----------
3  | TRIAL     | character | 0 (0.0%) |  HighTemperature | 60 (33.3%)
   |           |           |          |          Hypoxia | 60 (33.3%)
   |           |           |          |      LowSalinity | 60 (33.3%)
---+-----------+-----------+----------+------------------+-----------
4  | SMR_contr | numeric   | 0 (0.0%) |      [3.98, 6.7] |        180
---+-----------+-----------+----------+------------------+-----------
5  | CHANGE    | numeric   | 0 (0.0%) | [-55.01, 116.33] |        180
---------------------------------------------------------------------

norin |> modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
FISHID	60	0	30.5	17.4	1.0	30.5	60.0
MASS	175	0	40.2	8.0	22.5	39.3	64.1
SMR_contr	60	0	5.2	0.6	4.0	5.1	6.7
CHANGE	180	0	19.6	33.8	-55.0	16.4	116.3
TRIAL	N	%
HighTemperature	60	33.3
Hypoxia	60	33.3
LowSalinity	60	33.3

norin |> modelsummary::datasummary_skim(by = "TRIAL")

	TRIAL	Unique	Missing Pct.	Mean	SD	Min	Median	Max
FISHID	LowSalinity	60	0	30.5	17.5	1.0	30.5	60.0
	HighTemperature	60	0	30.5	17.5	1.0	30.5	60.0
	Hypoxia	60	0	30.5	17.5	1.0	30.5	60.0
MASS	LowSalinity	60	0	38.1	7.2	23.8	37.6	53.0
	HighTemperature	60	0	38.3	7.2	22.5	37.7	53.6
	Hypoxia	60	0	44.2	8.1	28.5	44.0	64.1
SMR_contr	LowSalinity	60	0	5.2	0.6	4.0	5.1	6.7
	HighTemperature	60	0	5.2	0.6	4.0	5.1	6.7
	Hypoxia	60	0	5.2	0.6	4.0	5.1	6.7
CHANGE	LowSalinity	60	0	8.2	34.4	-55.0	6.4	99.3
	HighTemperature	60	0	50.4	22.8	7.0	52.2	116.3
	Hypoxia	60	0	0.2	17.1	-41.2	-2.7	46.6
TRIAL	N	%
HighTemperature	60	33.3
Hypoxia	60	33.3
LowSalinity	60	33.3

4 Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}\\ \beta_0 \sim{} \mathcal{N}(16, 35)\\ \beta_{1-6} \sim{} \mathcal{N}(0, 70)\\ \]

where $\boldsymbol{\beta}$ and $\boldsymbol{\gamma}$ are vectors of the fixed and random effects parameters respectively and $\bf{X}$ is the model matrix representing the overall intercept and effects of temperature and (centered) mean fish size on SDA peak. $\bf{Z}$ represents a cell means model matrix for the random intercepts associated with individual fish.

Each of the fish were exposed to each of the three trials (treatments). So lets explore the distributions of responses within each of these trials.

ggplot(norin, aes(y = CHANGE, x = TRIAL)) +
  geom_boxplot()

Conclusions:

these seem reasonable enough

The researchers considered that physiological plasticity might be effected by the fishes basal metabolic rate. For example, a fish with relatively high metabolism might have less scope to increase this metabolism than a fish with a lower metabolism. Therefore, the SMR under control conditions (just prior to the specific trial) could be added as a continuous covariate.

Doing so would introduce two additional assumptions:

linearity: that the relationship between the response and control SMR is linear (if we intend to model it as a linear trend)
homogeneity of slopes: that the rate of change in response associated with a one unit change in control SMR is similar for each trial (unless we include an interaction)

ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
  geom_smooth(method = "lm") +
  geom_point()

`geom_smooth()` using formula = 'y ~ x'

ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
  geom_smooth() +
  geom_point()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

It might also be worth exploring how consistent the trial effect is across the different fish. This can give us an idea of whether the addition of a random intercept/slope model would be useful.

ggplot(norin, aes(y = CHANGE, x = as.numeric(FISHID), color = TRIAL)) +
  geom_point() +
  geom_line()

Conclusions:

generally, the difference between control and trial SMR is greatest for the High temperature trial.
there is less consistency about the relative effects of Hypoxia and Low Salinity.
hence a random intercept/slope model might be useful.

Finally, as an acknowledgement that the metabolic response might be influenced

by the mass of the fish, the researchers contemplated including fish Mass as a continuous covariate. There are numerous alternative ways to incorporate covariates that are known to impact on a response.

add them as an offset. This will include the covariate in the model, yet will not estimate a partial slope for this predictor. Rather the slope is assumed to be 1. Since the parameter does not need to be estimated, adding a predictor as an offset does not consume any degrees of freedom. However, it does assume that there is a slope of one. This approach can be useful as an alternative to using the response divided by the covariate as the response. That is, it is equivalent to using a response that has been standardised by the covariate, yet still allowing the modelling distribution that would be appropriate for the original response.
add them as a regular covariate. This will include the covariate in the model and estimate a partial slope parameter just like any other term. Although it will incur a degrees of freedom penalty, it does not assume a slope of 1 and could be modelled as non-linear if appropriate.

Lets explore the relationship between the response and Mass separately for each Trial.

ggplot(norin, aes(y = CHANGE, x = MASS, color = TRIAL)) +
  geom_point() +
  geom_smooth(method = "lm")

`geom_smooth()` using formula = 'y ~ x'

Conclusions:

there does not seem to be much of a relationship between the response and Mass. Nevertheless, in the presence of control SMN, it might still be a useful covariate at explaining some of the unexplained variability (and thus increasing the power of the model).
there is no evidence that the response and Mass are related in a non-linear manner.

5 Fit the model

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

norin.rstanarm <- stan_glmer(CHANGE ~ (1 | FISHID) + TRIAL * scale(SMR_contr, scale = FALSE) + scale(MASS, scale = FALSE),
  data = norin,
  family = gaussian(),
  iter = 5000,
  warmup = 2000,
  chains = 3,
  thin = 5,
  refresh = 0
)

norin.rstanarm |> prior_summary()

Priors for model 'norin.rstanarm' 
------
Intercept (after predictors centered)
  Specified prior:
    ~ normal(location = 20, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 20, scale = 85)

Coefficients
  Specified prior:
    ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
  Adjusted prior:
    ~ normal(location = [0,0,0,...], scale = [178.99,178.99,135.15,...])

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.03)

Covariance
 ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details

This tells us:

for the intercept (when the family is Gaussian), a normal prior with a mean of 20 and a standard deviation of 85 is used. A value of 2.5 is used for scaling all parameter standard deviations. The value of 20 is based on the mean of the response (mean(norin$CHANGE)) and the scaled standard deviation of 85 is based on multiplying the scaling factor by the standard deviation of the response.

2.5 * sd(norin$CHANGE)

[1] 84.6135

for the coefficients (in this case, just the difference between strong and weak innoculation), the default prior is a normal prior centered around 0 with a standard deviations of 178.99, 178.99, 135.15, 10.60, 34.2 and 34.2 repectively. This is then adjusted for the scale of the data by multiplying the 2.5 by the ratio of the standard deviation of the response by the stanard devation of the numerical dummy variables for the predictor (then rounded).

2.5 * sd(norin$CHANGE) / apply(model.matrix(~ TRIAL * SMR_contr + MASS, norin), 2, sd)

               (Intercept)               TRIALHypoxia 
                       Inf                  178.99305 
          TRIALLowSalinity                  SMR_contr 
                 178.99305                  135.14517 
                      MASS     TRIALHypoxia:SMR_contr 
                  10.59372                   34.19156 
TRIALLowSalinity:SMR_contr 
                  34.19156

the auxillary prior represents whatever additional prior is required for the nominated model. In the case of a Gaussian model, this will be $\sigma$, for negative binomial, it will be the reciprocal of dispersion, for gamma, it will be shape, etc . By default in rstanarm, this is a exponential with a rate of 1 which is then adjusted by devision with the standard deviation of the response.

1 / sd(norin$CHANGE)

[1] 0.02954611

Lets now run with priors only so that we can explore the range of values they allow in the posteriors.

norin.rstanarm1 <- update(norin.rstanarm, prior_PD = TRUE)

norin.rstanarm1 |>
  ggpredict(~ SMR_contr * TRIAL) |>
  plot(show_data = TRUE)

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

# OR
norin.rstanarm1 |>
  ggemmeans(~ SMR_contr * TRIAL) |>
  plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

Conclusions:

we see that the range of predictions is faily wide and the predicted means could range from a small negative number to a relatively large positive number.

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

$\beta_0$: normal centered at 17 with a standard deviation of 35
- mean of 17: since median(norin$CHANGE)
- sd of 35: since mad(norin$CHANGE)
$\beta_1$: normal centred at 0 with a standard deviation of 70
- sd of 70: since sd(norin$CHANGE) / apply(model.matrix(~scale(SMR_contr)*TRIAL + scale(MASS), norin), 2, sd)
$\sigma$: exponential with rate of 0.03
- since 1 / sd(norin$CHANGE)
$\Sigma$: decov with:
- regularization: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
- concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
- shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

I will also overlay the raw data for comparison.

norin.rstanarm2 <- stan_glmer(CHANGE ~ (1 | FISHID) + TRIAL * scale(SMR_contr, scale = FALSE) + offset(MASS),
  data = norin,
  family = gaussian(),
  prior_intercept = normal(17, 35, autoscale = FALSE),
  prior = normal(0, 70, autoscale = FALSE),
  prior_aux = rstanarm::exponential(0.03, autoscale = FALSE),
  prior_covariance = decov(1, 1, 1, 1),
  prior_PD = TRUE,
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0
)

norin.rstanarm2 |>
  ggpredict(~ SMR_contr * TRIAL) |>
  plot(show_data = TRUE)

Warning in sweep(eta, 2L, offset, `+`): STATS is longer than the extent of
'dim(x)[MARGIN]'

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

Now lets refit, conditioning on the data.

norin.rstanarm3 <- update(norin.rstanarm2, prior_PD = FALSE)

norin.rstanarm4 <- stan_glmer(CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + offset(MASS),
  data = norin,
  family = gaussian(),
  prior_intercept = normal(17, 35, autoscale = FALSE),
  prior = normal(0, 70, autoscale = FALSE),
  prior_aux = rstanarm::exponential(0.03, autoscale = FALSE),
  prior_covariance = decov(1, 1, 1, 1),
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0
)

Warning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
https://mc-stan.org/misc/warnings.html#bfmi-low

Warning: Examine the pairs() plot to diagnose sampling problems

preds <- norin.rstanarm4 |> posterior_predict(ndraws = 250, summary = FALSE)
norin.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = norin$CHANGE,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = FALSE
)
plot(norin.resids, quantreg = FALSE)

## Clearly an issue here!

get_variables(norin.rstanarm3)

 [1] "(Intercept)"                                     
 [2] "TRIALHypoxia"                                    
 [3] "TRIALLowSalinity"                                
 [4] "scale(SMR_contr, scale = FALSE)"                 
 [5] "TRIALHypoxia:scale(SMR_contr, scale = FALSE)"    
 [6] "TRIALLowSalinity:scale(SMR_contr, scale = FALSE)"
 [7] "b[(Intercept) FISHID:1]"                         
 [8] "b[(Intercept) FISHID:2]"                         
 [9] "b[(Intercept) FISHID:3]"                         
[10] "b[(Intercept) FISHID:4]"                         
[11] "b[(Intercept) FISHID:5]"                         
[12] "b[(Intercept) FISHID:6]"                         
[13] "b[(Intercept) FISHID:7]"                         
[14] "b[(Intercept) FISHID:8]"                         
[15] "b[(Intercept) FISHID:9]"                         
[16] "b[(Intercept) FISHID:10]"                        
[17] "b[(Intercept) FISHID:11]"                        
[18] "b[(Intercept) FISHID:12]"                        
[19] "b[(Intercept) FISHID:13]"                        
[20] "b[(Intercept) FISHID:14]"                        
[21] "b[(Intercept) FISHID:15]"                        
[22] "b[(Intercept) FISHID:16]"                        
[23] "b[(Intercept) FISHID:17]"                        
[24] "b[(Intercept) FISHID:18]"                        
[25] "b[(Intercept) FISHID:19]"                        
[26] "b[(Intercept) FISHID:20]"                        
[27] "b[(Intercept) FISHID:21]"                        
[28] "b[(Intercept) FISHID:22]"                        
[29] "b[(Intercept) FISHID:23]"                        
[30] "b[(Intercept) FISHID:24]"                        
[31] "b[(Intercept) FISHID:25]"                        
[32] "b[(Intercept) FISHID:26]"                        
[33] "b[(Intercept) FISHID:27]"                        
[34] "b[(Intercept) FISHID:28]"                        
[35] "b[(Intercept) FISHID:29]"                        
[36] "b[(Intercept) FISHID:30]"                        
[37] "b[(Intercept) FISHID:31]"                        
[38] "b[(Intercept) FISHID:32]"                        
[39] "b[(Intercept) FISHID:33]"                        
[40] "b[(Intercept) FISHID:34]"                        
[41] "b[(Intercept) FISHID:35]"                        
[42] "b[(Intercept) FISHID:36]"                        
[43] "b[(Intercept) FISHID:37]"                        
[44] "b[(Intercept) FISHID:38]"                        
[45] "b[(Intercept) FISHID:39]"                        
[46] "b[(Intercept) FISHID:40]"                        
[47] "b[(Intercept) FISHID:41]"                        
[48] "b[(Intercept) FISHID:42]"                        
[49] "b[(Intercept) FISHID:43]"                        
[50] "b[(Intercept) FISHID:44]"                        
[51] "b[(Intercept) FISHID:45]"                        
[52] "b[(Intercept) FISHID:46]"                        
[53] "b[(Intercept) FISHID:47]"                        
[54] "b[(Intercept) FISHID:48]"                        
[55] "b[(Intercept) FISHID:49]"                        
[56] "b[(Intercept) FISHID:50]"                        
[57] "b[(Intercept) FISHID:51]"                        
[58] "b[(Intercept) FISHID:52]"                        
[59] "b[(Intercept) FISHID:53]"                        
[60] "b[(Intercept) FISHID:54]"                        
[61] "b[(Intercept) FISHID:55]"                        
[62] "b[(Intercept) FISHID:56]"                        
[63] "b[(Intercept) FISHID:57]"                        
[64] "b[(Intercept) FISHID:58]"                        
[65] "b[(Intercept) FISHID:59]"                        
[66] "b[(Intercept) FISHID:60]"                        
[67] "sigma"                                           
[68] "Sigma[FISHID:(Intercept),(Intercept)]"           
[69] "accept_stat__"                                   
[70] "stepsize__"                                      
[71] "treedepth__"                                     
[72] "n_leapfrog__"                                    
[73] "divergent__"                                     
[74] "energy__"

posterior_vs_prior(norin.rstanarm3,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y"),
  regex_pars = "^.Intercept|TRIAL|SMR_contr|MASS|sigma|Sigma"
)


Drawing from prior...

Conclusions:

in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

norin.rstanarm3 |>
  ggpredict(~ SMR_contr * TRIAL) |>
  plot(show_data = TRUE)

Warning in sweep(eta, 2L, offset, `+`): STATS is longer than the extent of
'dim(x)[MARGIN]'

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

## OR
norin.rstanarm3 |>
  ggemmeans(~ SMR_contr * TRIAL) |>
  plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

norin.form <- bf(CHANGE ~  (1|FISHID)+TRIAL*SMR_contr+offset(MASS),
                   family = gaussian() 
                   )
options(width=100)
norin.form |> get_prior(data=norin)

                                 prior     class                       coef  group resp dpar nlpar
                                (flat)         b                                                  
                                (flat)         b                  SMR_contr                       
                                (flat)         b               TRIALHypoxia                       
                                (flat)         b     TRIALHypoxia:SMR_contr                       
                                (flat)         b           TRIALLowSalinity                       
                                (flat)         b TRIALLowSalinity:SMR_contr                       
 student_t(3, -23.7797222222222, 34.5) Intercept                                                  
                 student_t(3, 0, 34.5)        sd                                                  
                 student_t(3, 0, 34.5)        sd                            FISHID                
                 student_t(3, 0, 34.5)        sd                  Intercept FISHID                
                 student_t(3, 0, 34.5)     sigma                                                  
 lb ub       source
            default
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
            default
  0         default
  0    (vectorized)
  0    (vectorized)
  0         default

norin.brm <- brm(norin.form,
                 data=norin,
                 iter = 5000,
                 warmup = 1000,
                 chains = 3, cores = 3,
                 thin = 5,
                 refresh = 0,
                 backend = "rstan")

Compiling Stan program...

Start sampling

This tells us:

for the intercept (when the family is Gaussian), a t-distribution prior with a mean of 16.4 and a standard deviation of 34.5 is used:
- mean of 16.4, since median(norin$CHANGE)
- sd of 34.5 since mad(norin$CHANGE)
brms uses flat (inproper) priors for all population effects
the prior on sigma is a half t-distribution with a mean of 0 and standard deviation of 34.5:
- sd of 34.5 since mad(norin$CHANGE)
the hyperprior on the standard deviation is a half t-distribution with a mean of 0 and a standard deviation of 34.5.

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

$\beta_0$: normal centred at 16 with a standard deviation of 35
- mean of 16: since median(norin$CHANGE)
- sd of 35: since mad(norin$CHANGE)
$\beta_{1-2}$ (TRIAL): normal centred at 0 with a standard deviation of 70
- sd of 70: since sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
$\beta_{3}$ (SMR_contr): normal centred at 0 with a standard deviation of 54
- sd of 54: since sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
$\beta_{4}$ (MASS): normal centred at 0 with a standard deviation of 4
- sd of 4: since sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
$\beta_{5-6}$ (interaction): normal centred at 0 with a standard deviation of 14
- sd of 14: since sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
$\sigma$: half-cauchy with parameters 0 and 6 or else a student_t with parameters 0 and 34
- since sd(norin$CHANGE)
$\sigma_j$: half-cauchy with parameters 0 and 5.8
- since sqrt(sd(norin$CHANGE))
- we want this prior to have most mass close to zero for the purpose of regularisation
$\Sigma$: decov with:
- regularization: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
- concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
- shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

Note, for hierarchical models, the model will tend to want to have a large $sigma$ in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:

half-t: as the degrees of freedom approach infinity, this will approach a half-normal
half-cauchy: this is essentially a half-t with 1 degree of freedom
exponential

I will also overlay the raw data for comparison.

norin |>
  group_by(TRIAL) |>
  summarise(
    median(CHANGE),
    mad(CHANGE)
  )

# A tibble: 3 × 3
  TRIAL           `median(CHANGE)` `mad(CHANGE)`
  <fct>                      <dbl>         <dbl>
1 HighTemperature            52.2           23.2
2 Hypoxia                    -2.66          13.7
3 LowSalinity                 6.43          38.4

norin |>
  group_by(TRIAL, FISHID) |>
  summarise(
    median = median(CHANGE),
    MAD = mad(CHANGE)
  ) |>
  ungroup(FISHID) |>
  summarise(sd(median))

`summarise()` has grouped output by 'TRIAL'. You can override using the
`.groups` argument.

# A tibble: 3 × 2
  TRIAL           `sd(median)`
  <fct>                  <dbl>
1 HighTemperature         22.8
2 Hypoxia                 17.1
3 LowSalinity             34.4

sd(norin$CHANGE) / apply(model.matrix(~ TRIAL * SMR_contr + MASS, norin)[, -1], 2, sd)

              TRIALHypoxia           TRIALLowSalinity 
                 71.597218                  71.597218 
                 SMR_contr                       MASS 
                 54.058070                   4.237486 
    TRIALHypoxia:SMR_contr TRIALLowSalinity:SMR_contr 
                 13.676626                  13.676626

standist::visualize("normal(16,35)", xlim = c(-10, 100))

standist::visualize("normal(0, 70)", "normal(0, 54)", xlim = c(-200, 200))

standist::visualize("gamma(2, 1)", "gamma(35, 1)",
  "student_t(3,0, 34)",
  "cauchy(0, 5.8)",
  xlim = c(-10, 100)
)

priors <- prior(normal(16, 35), class = "Intercept") +
  prior(normal(0, 70), class = "b", coef = "TRIALHypoxia") +
  prior(normal(0, 70), class = "b", coef = "TRIALLowSalinity") +
  prior(normal(0, 54), class = "b", coef = "SMR_contr") +
  prior(normal(0, 4), class = "b", coef = "MASS") +
  prior(normal(0, 13), class = "b") +
  prior(student_t(3, 0, 34), class = "sigma") +
  ## prior(gamma(35, 1), class = 'sigma') +
  prior(cauchy(0, 5.8), class = "sd")
norin.form <- bf(CHANGE ~ (1 | FISHID) + TRIAL * SMR_contr + MASS,
  family = gaussian()
)
norin.brm2 <- brm(norin.form,
  data = norin,
  prior = priors,
  sample_prior = "only",
  iter = 5000,
  warmup = 1000,
  chains = 3, cores = 3,
  thin = 5,
  refresh = 0,
  backend = "rstan"
)

Compiling Stan program...

Start sampling

Note in the above model, the output may have included a warning message alerting us the presence of divergent transitions. Divergent transitions are an indication that the sampler has encountered poor sampling conditions - the more divergent transitions, the more severe the issue.

Typically, divergent transitions are the result of either:

a miss-specified model
priors that permit the sampler to drift into unsupported areas
complex posterior “features” (with high degrees of curvature) for which the sampler was inadequately tuned during the warmup phase

Accordingly, these divergent transitions can be addressed by either:

reviewing the model structure
adopting tighter priors
increase the adaptive delta from the default of 0.8 to closer to 1. The adaptive delta defines the average acceptance probability that the sampler should aspire to during the warmup phase. Increasing the adaptive delta results in a smaller step size (and thus fewer divergences and more robust samples) however it will also result in slower sampling speeds.

norin.brm2 |>
  ggpredict(~ SMR_contr * TRIAL) |>
  plot(show_data = TRUE)

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

norin.brm3 <- update(norin.brm2,
  sample_prior = "yes",
  control = list(adapt_delta = 0.99),
  refresh = 0,
  cores = 3
)

The desired updates require recompiling the model

Compiling Stan program...

Start sampling

norin.brm3 |>
  ggpredict(~ SMR_contr * TRIAL) |>
  plot(show_data = TRUE)

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

norin.brm3 |> get_variables()

 [1] "b_Intercept"                        "b_TRIALHypoxia"                    
 [3] "b_TRIALLowSalinity"                 "b_SMR_contr"                       
 [5] "b_MASS"                             "b_TRIALHypoxia:SMR_contr"          
 [7] "b_TRIALLowSalinity:SMR_contr"       "sd_FISHID__Intercept"              
 [9] "sigma"                              "Intercept"                         
[11] "r_FISHID[1,Intercept]"              "r_FISHID[2,Intercept]"             
[13] "r_FISHID[3,Intercept]"              "r_FISHID[4,Intercept]"             
[15] "r_FISHID[5,Intercept]"              "r_FISHID[6,Intercept]"             
[17] "r_FISHID[7,Intercept]"              "r_FISHID[8,Intercept]"             
[19] "r_FISHID[9,Intercept]"              "r_FISHID[10,Intercept]"            
[21] "r_FISHID[11,Intercept]"             "r_FISHID[12,Intercept]"            
[23] "r_FISHID[13,Intercept]"             "r_FISHID[14,Intercept]"            
[25] "r_FISHID[15,Intercept]"             "r_FISHID[16,Intercept]"            
[27] "r_FISHID[17,Intercept]"             "r_FISHID[18,Intercept]"            
[29] "r_FISHID[19,Intercept]"             "r_FISHID[20,Intercept]"            
[31] "r_FISHID[21,Intercept]"             "r_FISHID[22,Intercept]"            
[33] "r_FISHID[23,Intercept]"             "r_FISHID[24,Intercept]"            
[35] "r_FISHID[25,Intercept]"             "r_FISHID[26,Intercept]"            
[37] "r_FISHID[27,Intercept]"             "r_FISHID[28,Intercept]"            
[39] "r_FISHID[29,Intercept]"             "r_FISHID[30,Intercept]"            
[41] "r_FISHID[31,Intercept]"             "r_FISHID[32,Intercept]"            
[43] "r_FISHID[33,Intercept]"             "r_FISHID[34,Intercept]"            
[45] "r_FISHID[35,Intercept]"             "r_FISHID[36,Intercept]"            
[47] "r_FISHID[37,Intercept]"             "r_FISHID[38,Intercept]"            
[49] "r_FISHID[39,Intercept]"             "r_FISHID[40,Intercept]"            
[51] "r_FISHID[41,Intercept]"             "r_FISHID[42,Intercept]"            
[53] "r_FISHID[43,Intercept]"             "r_FISHID[44,Intercept]"            
[55] "r_FISHID[45,Intercept]"             "r_FISHID[46,Intercept]"            
[57] "r_FISHID[47,Intercept]"             "r_FISHID[48,Intercept]"            
[59] "r_FISHID[49,Intercept]"             "r_FISHID[50,Intercept]"            
[61] "r_FISHID[51,Intercept]"             "r_FISHID[52,Intercept]"            
[63] "r_FISHID[53,Intercept]"             "r_FISHID[54,Intercept]"            
[65] "r_FISHID[55,Intercept]"             "r_FISHID[56,Intercept]"            
[67] "r_FISHID[57,Intercept]"             "r_FISHID[58,Intercept]"            
[69] "r_FISHID[59,Intercept]"             "r_FISHID[60,Intercept]"            
[71] "prior_Intercept"                    "prior_b_TRIALHypoxia"              
[73] "prior_b_TRIALLowSalinity"           "prior_b_SMR_contr"                 
[75] "prior_b_MASS"                       "prior_b_TRIALHypoxia:SMR_contr"    
[77] "prior_b_TRIALLowSalinity:SMR_contr" "prior_sigma"                       
[79] "prior_sd_FISHID"                    "lprior"                            
[81] "lp__"                               "accept_stat__"                     
[83] "stepsize__"                         "treedepth__"                       
[85] "n_leapfrog__"                       "divergent__"                       
[87] "energy__"

norin.brm3 |>
  hypothesis("TRIALHypoxia=0") |>
  plot()

norin.brm3 |>
  hypothesis("SMR_contr=0") |>
  plot()

## norin.brm3 |> SUYR_prior_and_posterior()

While we are here, we might like to explore a random intercept/slope model

priors <- prior(normal(16, 35), class = "Intercept") +
  prior(normal(0, 70), class = "b", coef = "TRIALHypoxia") +
  prior(normal(0, 70), class = "b", coef = "TRIALLowSalinity") +
  prior(normal(0, 54), class = "b", coef = "SMR_contr") +
  prior(normal(0, 4), class = "b", coef = "MASS") +
  ## prior(normal(0, 70), class = 'b') +
  prior(student_t(3, 0, 34), class = "sigma") +
  prior(cauchy(0, 5.8), class = "sd") +
  prior(lkj_corr_cholesky(1), class = "cor")
norin.form <- bf(CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + MASS,
  ## sigma ~ TRIAL*SMR_contr + offset(MASS) + (1|FISHID),
  family = gaussian()
)
norin.brm4 <- brm(norin.form,
  data = norin,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 1000,
  chains = 3, cores = 3,
  thin = 5,
  refresh = 0,
  control = list(adapt_delta = 0.99),
  backend = "rstan"
)

Compiling Stan program...

Start sampling

Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
https://mc-stan.org/misc/warnings.html#bfmi-low

Warning: Examine the pairs() plot to diagnose sampling problems

norin.brm4 |> get_variables()

  [1] "b_Intercept"                               
  [2] "b_TRIALHypoxia"                            
  [3] "b_TRIALLowSalinity"                        
  [4] "b_SMR_contr"                               
  [5] "b_MASS"                                    
  [6] "b_TRIALHypoxia:SMR_contr"                  
  [7] "b_TRIALLowSalinity:SMR_contr"              
  [8] "sd_FISHID__Intercept"                      
  [9] "sd_FISHID__TRIALHypoxia"                   
 [10] "sd_FISHID__TRIALLowSalinity"               
 [11] "cor_FISHID__Intercept__TRIALHypoxia"       
 [12] "cor_FISHID__Intercept__TRIALLowSalinity"   
 [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
 [14] "sigma"                                     
 [15] "Intercept"                                 
 [16] "r_FISHID[1,Intercept]"                     
 [17] "r_FISHID[2,Intercept]"                     
 [18] "r_FISHID[3,Intercept]"                     
 [19] "r_FISHID[4,Intercept]"                     
 [20] "r_FISHID[5,Intercept]"                     
 [21] "r_FISHID[6,Intercept]"                     
 [22] "r_FISHID[7,Intercept]"                     
 [23] "r_FISHID[8,Intercept]"                     
 [24] "r_FISHID[9,Intercept]"                     
 [25] "r_FISHID[10,Intercept]"                    
 [26] "r_FISHID[11,Intercept]"                    
 [27] "r_FISHID[12,Intercept]"                    
 [28] "r_FISHID[13,Intercept]"                    
 [29] "r_FISHID[14,Intercept]"                    
 [30] "r_FISHID[15,Intercept]"                    
 [31] "r_FISHID[16,Intercept]"                    
 [32] "r_FISHID[17,Intercept]"                    
 [33] "r_FISHID[18,Intercept]"                    
 [34] "r_FISHID[19,Intercept]"                    
 [35] "r_FISHID[20,Intercept]"                    
 [36] "r_FISHID[21,Intercept]"                    
 [37] "r_FISHID[22,Intercept]"                    
 [38] "r_FISHID[23,Intercept]"                    
 [39] "r_FISHID[24,Intercept]"                    
 [40] "r_FISHID[25,Intercept]"                    
 [41] "r_FISHID[26,Intercept]"                    
 [42] "r_FISHID[27,Intercept]"                    
 [43] "r_FISHID[28,Intercept]"                    
 [44] "r_FISHID[29,Intercept]"                    
 [45] "r_FISHID[30,Intercept]"                    
 [46] "r_FISHID[31,Intercept]"                    
 [47] "r_FISHID[32,Intercept]"                    
 [48] "r_FISHID[33,Intercept]"                    
 [49] "r_FISHID[34,Intercept]"                    
 [50] "r_FISHID[35,Intercept]"                    
 [51] "r_FISHID[36,Intercept]"                    
 [52] "r_FISHID[37,Intercept]"                    
 [53] "r_FISHID[38,Intercept]"                    
 [54] "r_FISHID[39,Intercept]"                    
 [55] "r_FISHID[40,Intercept]"                    
 [56] "r_FISHID[41,Intercept]"                    
 [57] "r_FISHID[42,Intercept]"                    
 [58] "r_FISHID[43,Intercept]"                    
 [59] "r_FISHID[44,Intercept]"                    
 [60] "r_FISHID[45,Intercept]"                    
 [61] "r_FISHID[46,Intercept]"                    
 [62] "r_FISHID[47,Intercept]"                    
 [63] "r_FISHID[48,Intercept]"                    
 [64] "r_FISHID[49,Intercept]"                    
 [65] "r_FISHID[50,Intercept]"                    
 [66] "r_FISHID[51,Intercept]"                    
 [67] "r_FISHID[52,Intercept]"                    
 [68] "r_FISHID[53,Intercept]"                    
 [69] "r_FISHID[54,Intercept]"                    
 [70] "r_FISHID[55,Intercept]"                    
 [71] "r_FISHID[56,Intercept]"                    
 [72] "r_FISHID[57,Intercept]"                    
 [73] "r_FISHID[58,Intercept]"                    
 [74] "r_FISHID[59,Intercept]"                    
 [75] "r_FISHID[60,Intercept]"                    
 [76] "r_FISHID[1,TRIALHypoxia]"                  
 [77] "r_FISHID[2,TRIALHypoxia]"                  
 [78] "r_FISHID[3,TRIALHypoxia]"                  
 [79] "r_FISHID[4,TRIALHypoxia]"                  
 [80] "r_FISHID[5,TRIALHypoxia]"                  
 [81] "r_FISHID[6,TRIALHypoxia]"                  
 [82] "r_FISHID[7,TRIALHypoxia]"                  
 [83] "r_FISHID[8,TRIALHypoxia]"                  
 [84] "r_FISHID[9,TRIALHypoxia]"                  
 [85] "r_FISHID[10,TRIALHypoxia]"                 
 [86] "r_FISHID[11,TRIALHypoxia]"                 
 [87] "r_FISHID[12,TRIALHypoxia]"                 
 [88] "r_FISHID[13,TRIALHypoxia]"                 
 [89] "r_FISHID[14,TRIALHypoxia]"                 
 [90] "r_FISHID[15,TRIALHypoxia]"                 
 [91] "r_FISHID[16,TRIALHypoxia]"                 
 [92] "r_FISHID[17,TRIALHypoxia]"                 
 [93] "r_FISHID[18,TRIALHypoxia]"                 
 [94] "r_FISHID[19,TRIALHypoxia]"                 
 [95] "r_FISHID[20,TRIALHypoxia]"                 
 [96] "r_FISHID[21,TRIALHypoxia]"                 
 [97] "r_FISHID[22,TRIALHypoxia]"                 
 [98] "r_FISHID[23,TRIALHypoxia]"                 
 [99] "r_FISHID[24,TRIALHypoxia]"                 
[100] "r_FISHID[25,TRIALHypoxia]"                 
[101] "r_FISHID[26,TRIALHypoxia]"                 
[102] "r_FISHID[27,TRIALHypoxia]"                 
[103] "r_FISHID[28,TRIALHypoxia]"                 
[104] "r_FISHID[29,TRIALHypoxia]"                 
[105] "r_FISHID[30,TRIALHypoxia]"                 
[106] "r_FISHID[31,TRIALHypoxia]"                 
[107] "r_FISHID[32,TRIALHypoxia]"                 
[108] "r_FISHID[33,TRIALHypoxia]"                 
[109] "r_FISHID[34,TRIALHypoxia]"                 
[110] "r_FISHID[35,TRIALHypoxia]"                 
[111] "r_FISHID[36,TRIALHypoxia]"                 
[112] "r_FISHID[37,TRIALHypoxia]"                 
[113] "r_FISHID[38,TRIALHypoxia]"                 
[114] "r_FISHID[39,TRIALHypoxia]"                 
[115] "r_FISHID[40,TRIALHypoxia]"                 
[116] "r_FISHID[41,TRIALHypoxia]"                 
[117] "r_FISHID[42,TRIALHypoxia]"                 
[118] "r_FISHID[43,TRIALHypoxia]"                 
[119] "r_FISHID[44,TRIALHypoxia]"                 
[120] "r_FISHID[45,TRIALHypoxia]"                 
[121] "r_FISHID[46,TRIALHypoxia]"                 
[122] "r_FISHID[47,TRIALHypoxia]"                 
[123] "r_FISHID[48,TRIALHypoxia]"                 
[124] "r_FISHID[49,TRIALHypoxia]"                 
[125] "r_FISHID[50,TRIALHypoxia]"                 
[126] "r_FISHID[51,TRIALHypoxia]"                 
[127] "r_FISHID[52,TRIALHypoxia]"                 
[128] "r_FISHID[53,TRIALHypoxia]"                 
[129] "r_FISHID[54,TRIALHypoxia]"                 
[130] "r_FISHID[55,TRIALHypoxia]"                 
[131] "r_FISHID[56,TRIALHypoxia]"                 
[132] "r_FISHID[57,TRIALHypoxia]"                 
[133] "r_FISHID[58,TRIALHypoxia]"                 
[134] "r_FISHID[59,TRIALHypoxia]"                 
[135] "r_FISHID[60,TRIALHypoxia]"                 
[136] "r_FISHID[1,TRIALLowSalinity]"              
[137] "r_FISHID[2,TRIALLowSalinity]"              
[138] "r_FISHID[3,TRIALLowSalinity]"              
[139] "r_FISHID[4,TRIALLowSalinity]"              
[140] "r_FISHID[5,TRIALLowSalinity]"              
[141] "r_FISHID[6,TRIALLowSalinity]"              
[142] "r_FISHID[7,TRIALLowSalinity]"              
[143] "r_FISHID[8,TRIALLowSalinity]"              
[144] "r_FISHID[9,TRIALLowSalinity]"              
[145] "r_FISHID[10,TRIALLowSalinity]"             
[146] "r_FISHID[11,TRIALLowSalinity]"             
[147] "r_FISHID[12,TRIALLowSalinity]"             
[148] "r_FISHID[13,TRIALLowSalinity]"             
[149] "r_FISHID[14,TRIALLowSalinity]"             
[150] "r_FISHID[15,TRIALLowSalinity]"             
[151] "r_FISHID[16,TRIALLowSalinity]"             
[152] "r_FISHID[17,TRIALLowSalinity]"             
[153] "r_FISHID[18,TRIALLowSalinity]"             
[154] "r_FISHID[19,TRIALLowSalinity]"             
[155] "r_FISHID[20,TRIALLowSalinity]"             
[156] "r_FISHID[21,TRIALLowSalinity]"             
[157] "r_FISHID[22,TRIALLowSalinity]"             
[158] "r_FISHID[23,TRIALLowSalinity]"             
[159] "r_FISHID[24,TRIALLowSalinity]"             
[160] "r_FISHID[25,TRIALLowSalinity]"             
[161] "r_FISHID[26,TRIALLowSalinity]"             
[162] "r_FISHID[27,TRIALLowSalinity]"             
[163] "r_FISHID[28,TRIALLowSalinity]"             
[164] "r_FISHID[29,TRIALLowSalinity]"             
[165] "r_FISHID[30,TRIALLowSalinity]"             
[166] "r_FISHID[31,TRIALLowSalinity]"             
[167] "r_FISHID[32,TRIALLowSalinity]"             
[168] "r_FISHID[33,TRIALLowSalinity]"             
[169] "r_FISHID[34,TRIALLowSalinity]"             
[170] "r_FISHID[35,TRIALLowSalinity]"             
[171] "r_FISHID[36,TRIALLowSalinity]"             
[172] "r_FISHID[37,TRIALLowSalinity]"             
[173] "r_FISHID[38,TRIALLowSalinity]"             
[174] "r_FISHID[39,TRIALLowSalinity]"             
[175] "r_FISHID[40,TRIALLowSalinity]"             
[176] "r_FISHID[41,TRIALLowSalinity]"             
[177] "r_FISHID[42,TRIALLowSalinity]"             
[178] "r_FISHID[43,TRIALLowSalinity]"             
[179] "r_FISHID[44,TRIALLowSalinity]"             
[180] "r_FISHID[45,TRIALLowSalinity]"             
[181] "r_FISHID[46,TRIALLowSalinity]"             
[182] "r_FISHID[47,TRIALLowSalinity]"             
[183] "r_FISHID[48,TRIALLowSalinity]"             
[184] "r_FISHID[49,TRIALLowSalinity]"             
[185] "r_FISHID[50,TRIALLowSalinity]"             
[186] "r_FISHID[51,TRIALLowSalinity]"             
[187] "r_FISHID[52,TRIALLowSalinity]"             
[188] "r_FISHID[53,TRIALLowSalinity]"             
[189] "r_FISHID[54,TRIALLowSalinity]"             
[190] "r_FISHID[55,TRIALLowSalinity]"             
[191] "r_FISHID[56,TRIALLowSalinity]"             
[192] "r_FISHID[57,TRIALLowSalinity]"             
[193] "r_FISHID[58,TRIALLowSalinity]"             
[194] "r_FISHID[59,TRIALLowSalinity]"             
[195] "r_FISHID[60,TRIALLowSalinity]"             
[196] "prior_Intercept"                           
[197] "prior_b_TRIALHypoxia"                      
[198] "prior_b_TRIALLowSalinity"                  
[199] "prior_b_SMR_contr"                         
[200] "prior_b_MASS"                              
[201] "prior_sigma"                               
[202] "prior_sd_FISHID"                           
[203] "prior_cor_FISHID"                          
[204] "lprior"                                    
[205] "lp__"                                      
[206] "accept_stat__"                             
[207] "stepsize__"                                
[208] "treedepth__"                               
[209] "n_leapfrog__"                              
[210] "divergent__"                               
[211] "energy__"

norin.brm4 |>
  hypothesis("TRIALHypoxia=0") |>
  plot()

norin.brm4 |>
  hypothesis("SMR_contr=0") |>
  plot()

## norin.brm4 |> SUYR_prior_and_posterior()

norin.brm4 |>
  posterior_samples() |>
  dplyr::select(-`lp__`) |>
  pivot_longer(everything(), names_to = "key") |>
  filter(!str_detect(key, "^r")) |>
  mutate(
    Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
    Class = case_when(
      str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
      str_detect(key, "b_TRIAL.*|prior_b_TRIAL.*") & !str_detect(key, ".*\\:.*") ~ "TRIAL",
      str_detect(key, "b_SMR_contr|prior_b_SMR_contr") ~ "SMR",
      str_detect(key, "b_MASS|prior_b_MASS") ~ "MASS",
      str_detect(key, ".*\\:.*|prior_b_.*\\:.*") ~ "Interaction",
      str_detect(key, "sd") ~ "sd",
      str_detect(key, "^cor|prior_cor") ~ "cor",
      str_detect(key, "sigma") ~ "sigma"
    ),
    Par = str_replace(key, "b_", "")
  ) |>
  ggplot(aes(x = Type, y = value, color = Par)) +
  stat_pointinterval(position = position_dodge()) +
  facet_wrap(~Class, scales = "free")

Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.

We can compare the two models using LOO

(l.1 <- norin.brm3 |> loo())


Computed from 2400 by 180 log-likelihood matrix.

         Estimate   SE
elpd_loo   -832.1 13.7
p_loo        24.5  3.9
looic      1664.2 27.4
------
MCSE of elpd_loo is 0.2.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

(l.2 <- norin.brm4 |> loo())

Warning: Found 54 observations with a pareto_k > 0.7 in model 'norin.brm4'. We
recommend to set 'moment_match = TRUE' in order to perform moment matching for
problematic observations.


Computed from 2400 by 180 log-likelihood matrix.

         Estimate   SE
elpd_loo   -780.1  9.2
p_loo        85.7  6.1
looic      1560.3 18.4
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 0.7]).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     126   70.0%   47      
   (0.7, 1]   (bad)       47   26.1%   <NA>    
   (1, Inf)   (very bad)   7    3.9%   <NA>    
See help('pareto-k-diagnostic') for details.

loo_compare(l.1, l.2)

           elpd_diff se_diff
norin.brm4   0.0       0.0  
norin.brm3 -52.0      10.2

There is substantially more support for the more complex random intercept/slope model over the simpler random intercept only model. Consequently, we will continue with the random intercept/slope model.

6 MCMC sampling diagnostics

brms

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

See list of available diagnostics by name

available_mcmc()

bayesplot MCMC module:
  mcmc_acf
  mcmc_acf_bar
  mcmc_areas
  mcmc_areas_ridges
  mcmc_combo
  mcmc_dens
  mcmc_dens_chains
  mcmc_dens_overlay
  mcmc_hex
  mcmc_hist
  mcmc_hist_by_chain
  mcmc_intervals
  mcmc_neff
  mcmc_neff_hist
  mcmc_nuts_acceptance
  mcmc_nuts_divergence
  mcmc_nuts_energy
  mcmc_nuts_stepsize
  mcmc_nuts_treedepth
  mcmc_pairs
  mcmc_parcoord
  mcmc_rank_ecdf
  mcmc_rank_hist
  mcmc_rank_overlay
  mcmc_recover_hist
  mcmc_recover_intervals
  mcmc_recover_scatter
  mcmc_rhat
  mcmc_rhat_hist
  mcmc_scatter
  mcmc_trace
  mcmc_trace_highlight
  mcmc_violin

Of these, we will focus on:

trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

norin.brm4 |> mcmc_plot(type = "trace")

No divergences to plot.

The chains appear well mixed and very similar

acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

norin.brm4 |> mcmc_plot(type = "acf_bar")

There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!

rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

norin.brm4 |> mcmc_plot(type = "rhat_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

norin.brm4 |> mcmc_plot(type = "neff_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.

More diagnostics

norin.brm4 |> mcmc_plot(type = "combo")

norin.brm4 |> mcmc_plot(type = "violin")

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

norin.brm4 |> get_variables()

  [1] "b_Intercept"                               
  [2] "b_TRIALHypoxia"                            
  [3] "b_TRIALLowSalinity"                        
  [4] "b_SMR_contr"                               
  [5] "b_MASS"                                    
  [6] "b_TRIALHypoxia:SMR_contr"                  
  [7] "b_TRIALLowSalinity:SMR_contr"              
  [8] "sd_FISHID__Intercept"                      
  [9] "sd_FISHID__TRIALHypoxia"                   
 [10] "sd_FISHID__TRIALLowSalinity"               
 [11] "cor_FISHID__Intercept__TRIALHypoxia"       
 [12] "cor_FISHID__Intercept__TRIALLowSalinity"   
 [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
 [14] "sigma"                                     
 [15] "Intercept"                                 
 [16] "r_FISHID[1,Intercept]"                     
 [17] "r_FISHID[2,Intercept]"                     
 [18] "r_FISHID[3,Intercept]"                     
 [19] "r_FISHID[4,Intercept]"                     
 [20] "r_FISHID[5,Intercept]"                     
 [21] "r_FISHID[6,Intercept]"                     
 [22] "r_FISHID[7,Intercept]"                     
 [23] "r_FISHID[8,Intercept]"                     
 [24] "r_FISHID[9,Intercept]"                     
 [25] "r_FISHID[10,Intercept]"                    
 [26] "r_FISHID[11,Intercept]"                    
 [27] "r_FISHID[12,Intercept]"                    
 [28] "r_FISHID[13,Intercept]"                    
 [29] "r_FISHID[14,Intercept]"                    
 [30] "r_FISHID[15,Intercept]"                    
 [31] "r_FISHID[16,Intercept]"                    
 [32] "r_FISHID[17,Intercept]"                    
 [33] "r_FISHID[18,Intercept]"                    
 [34] "r_FISHID[19,Intercept]"                    
 [35] "r_FISHID[20,Intercept]"                    
 [36] "r_FISHID[21,Intercept]"                    
 [37] "r_FISHID[22,Intercept]"                    
 [38] "r_FISHID[23,Intercept]"                    
 [39] "r_FISHID[24,Intercept]"                    
 [40] "r_FISHID[25,Intercept]"                    
 [41] "r_FISHID[26,Intercept]"                    
 [42] "r_FISHID[27,Intercept]"                    
 [43] "r_FISHID[28,Intercept]"                    
 [44] "r_FISHID[29,Intercept]"                    
 [45] "r_FISHID[30,Intercept]"                    
 [46] "r_FISHID[31,Intercept]"                    
 [47] "r_FISHID[32,Intercept]"                    
 [48] "r_FISHID[33,Intercept]"                    
 [49] "r_FISHID[34,Intercept]"                    
 [50] "r_FISHID[35,Intercept]"                    
 [51] "r_FISHID[36,Intercept]"                    
 [52] "r_FISHID[37,Intercept]"                    
 [53] "r_FISHID[38,Intercept]"                    
 [54] "r_FISHID[39,Intercept]"                    
 [55] "r_FISHID[40,Intercept]"                    
 [56] "r_FISHID[41,Intercept]"                    
 [57] "r_FISHID[42,Intercept]"                    
 [58] "r_FISHID[43,Intercept]"                    
 [59] "r_FISHID[44,Intercept]"                    
 [60] "r_FISHID[45,Intercept]"                    
 [61] "r_FISHID[46,Intercept]"                    
 [62] "r_FISHID[47,Intercept]"                    
 [63] "r_FISHID[48,Intercept]"                    
 [64] "r_FISHID[49,Intercept]"                    
 [65] "r_FISHID[50,Intercept]"                    
 [66] "r_FISHID[51,Intercept]"                    
 [67] "r_FISHID[52,Intercept]"                    
 [68] "r_FISHID[53,Intercept]"                    
 [69] "r_FISHID[54,Intercept]"                    
 [70] "r_FISHID[55,Intercept]"                    
 [71] "r_FISHID[56,Intercept]"                    
 [72] "r_FISHID[57,Intercept]"                    
 [73] "r_FISHID[58,Intercept]"                    
 [74] "r_FISHID[59,Intercept]"                    
 [75] "r_FISHID[60,Intercept]"                    
 [76] "r_FISHID[1,TRIALHypoxia]"                  
 [77] "r_FISHID[2,TRIALHypoxia]"                  
 [78] "r_FISHID[3,TRIALHypoxia]"                  
 [79] "r_FISHID[4,TRIALHypoxia]"                  
 [80] "r_FISHID[5,TRIALHypoxia]"                  
 [81] "r_FISHID[6,TRIALHypoxia]"                  
 [82] "r_FISHID[7,TRIALHypoxia]"                  
 [83] "r_FISHID[8,TRIALHypoxia]"                  
 [84] "r_FISHID[9,TRIALHypoxia]"                  
 [85] "r_FISHID[10,TRIALHypoxia]"                 
 [86] "r_FISHID[11,TRIALHypoxia]"                 
 [87] "r_FISHID[12,TRIALHypoxia]"                 
 [88] "r_FISHID[13,TRIALHypoxia]"                 
 [89] "r_FISHID[14,TRIALHypoxia]"                 
 [90] "r_FISHID[15,TRIALHypoxia]"                 
 [91] "r_FISHID[16,TRIALHypoxia]"                 
 [92] "r_FISHID[17,TRIALHypoxia]"                 
 [93] "r_FISHID[18,TRIALHypoxia]"                 
 [94] "r_FISHID[19,TRIALHypoxia]"                 
 [95] "r_FISHID[20,TRIALHypoxia]"                 
 [96] "r_FISHID[21,TRIALHypoxia]"                 
 [97] "r_FISHID[22,TRIALHypoxia]"                 
 [98] "r_FISHID[23,TRIALHypoxia]"                 
 [99] "r_FISHID[24,TRIALHypoxia]"                 
[100] "r_FISHID[25,TRIALHypoxia]"                 
[101] "r_FISHID[26,TRIALHypoxia]"                 
[102] "r_FISHID[27,TRIALHypoxia]"                 
[103] "r_FISHID[28,TRIALHypoxia]"                 
[104] "r_FISHID[29,TRIALHypoxia]"                 
[105] "r_FISHID[30,TRIALHypoxia]"                 
[106] "r_FISHID[31,TRIALHypoxia]"                 
[107] "r_FISHID[32,TRIALHypoxia]"                 
[108] "r_FISHID[33,TRIALHypoxia]"                 
[109] "r_FISHID[34,TRIALHypoxia]"                 
[110] "r_FISHID[35,TRIALHypoxia]"                 
[111] "r_FISHID[36,TRIALHypoxia]"                 
[112] "r_FISHID[37,TRIALHypoxia]"                 
[113] "r_FISHID[38,TRIALHypoxia]"                 
[114] "r_FISHID[39,TRIALHypoxia]"                 
[115] "r_FISHID[40,TRIALHypoxia]"                 
[116] "r_FISHID[41,TRIALHypoxia]"                 
[117] "r_FISHID[42,TRIALHypoxia]"                 
[118] "r_FISHID[43,TRIALHypoxia]"                 
[119] "r_FISHID[44,TRIALHypoxia]"                 
[120] "r_FISHID[45,TRIALHypoxia]"                 
[121] "r_FISHID[46,TRIALHypoxia]"                 
[122] "r_FISHID[47,TRIALHypoxia]"                 
[123] "r_FISHID[48,TRIALHypoxia]"                 
[124] "r_FISHID[49,TRIALHypoxia]"                 
[125] "r_FISHID[50,TRIALHypoxia]"                 
[126] "r_FISHID[51,TRIALHypoxia]"                 
[127] "r_FISHID[52,TRIALHypoxia]"                 
[128] "r_FISHID[53,TRIALHypoxia]"                 
[129] "r_FISHID[54,TRIALHypoxia]"                 
[130] "r_FISHID[55,TRIALHypoxia]"                 
[131] "r_FISHID[56,TRIALHypoxia]"                 
[132] "r_FISHID[57,TRIALHypoxia]"                 
[133] "r_FISHID[58,TRIALHypoxia]"                 
[134] "r_FISHID[59,TRIALHypoxia]"                 
[135] "r_FISHID[60,TRIALHypoxia]"                 
[136] "r_FISHID[1,TRIALLowSalinity]"              
[137] "r_FISHID[2,TRIALLowSalinity]"              
[138] "r_FISHID[3,TRIALLowSalinity]"              
[139] "r_FISHID[4,TRIALLowSalinity]"              
[140] "r_FISHID[5,TRIALLowSalinity]"              
[141] "r_FISHID[6,TRIALLowSalinity]"              
[142] "r_FISHID[7,TRIALLowSalinity]"              
[143] "r_FISHID[8,TRIALLowSalinity]"              
[144] "r_FISHID[9,TRIALLowSalinity]"              
[145] "r_FISHID[10,TRIALLowSalinity]"             
[146] "r_FISHID[11,TRIALLowSalinity]"             
[147] "r_FISHID[12,TRIALLowSalinity]"             
[148] "r_FISHID[13,TRIALLowSalinity]"             
[149] "r_FISHID[14,TRIALLowSalinity]"             
[150] "r_FISHID[15,TRIALLowSalinity]"             
[151] "r_FISHID[16,TRIALLowSalinity]"             
[152] "r_FISHID[17,TRIALLowSalinity]"             
[153] "r_FISHID[18,TRIALLowSalinity]"             
[154] "r_FISHID[19,TRIALLowSalinity]"             
[155] "r_FISHID[20,TRIALLowSalinity]"             
[156] "r_FISHID[21,TRIALLowSalinity]"             
[157] "r_FISHID[22,TRIALLowSalinity]"             
[158] "r_FISHID[23,TRIALLowSalinity]"             
[159] "r_FISHID[24,TRIALLowSalinity]"             
[160] "r_FISHID[25,TRIALLowSalinity]"             
[161] "r_FISHID[26,TRIALLowSalinity]"             
[162] "r_FISHID[27,TRIALLowSalinity]"             
[163] "r_FISHID[28,TRIALLowSalinity]"             
[164] "r_FISHID[29,TRIALLowSalinity]"             
[165] "r_FISHID[30,TRIALLowSalinity]"             
[166] "r_FISHID[31,TRIALLowSalinity]"             
[167] "r_FISHID[32,TRIALLowSalinity]"             
[168] "r_FISHID[33,TRIALLowSalinity]"             
[169] "r_FISHID[34,TRIALLowSalinity]"             
[170] "r_FISHID[35,TRIALLowSalinity]"             
[171] "r_FISHID[36,TRIALLowSalinity]"             
[172] "r_FISHID[37,TRIALLowSalinity]"             
[173] "r_FISHID[38,TRIALLowSalinity]"             
[174] "r_FISHID[39,TRIALLowSalinity]"             
[175] "r_FISHID[40,TRIALLowSalinity]"             
[176] "r_FISHID[41,TRIALLowSalinity]"             
[177] "r_FISHID[42,TRIALLowSalinity]"             
[178] "r_FISHID[43,TRIALLowSalinity]"             
[179] "r_FISHID[44,TRIALLowSalinity]"             
[180] "r_FISHID[45,TRIALLowSalinity]"             
[181] "r_FISHID[46,TRIALLowSalinity]"             
[182] "r_FISHID[47,TRIALLowSalinity]"             
[183] "r_FISHID[48,TRIALLowSalinity]"             
[184] "r_FISHID[49,TRIALLowSalinity]"             
[185] "r_FISHID[50,TRIALLowSalinity]"             
[186] "r_FISHID[51,TRIALLowSalinity]"             
[187] "r_FISHID[52,TRIALLowSalinity]"             
[188] "r_FISHID[53,TRIALLowSalinity]"             
[189] "r_FISHID[54,TRIALLowSalinity]"             
[190] "r_FISHID[55,TRIALLowSalinity]"             
[191] "r_FISHID[56,TRIALLowSalinity]"             
[192] "r_FISHID[57,TRIALLowSalinity]"             
[193] "r_FISHID[58,TRIALLowSalinity]"             
[194] "r_FISHID[59,TRIALLowSalinity]"             
[195] "r_FISHID[60,TRIALLowSalinity]"             
[196] "prior_Intercept"                           
[197] "prior_b_TRIALHypoxia"                      
[198] "prior_b_TRIALLowSalinity"                  
[199] "prior_b_SMR_contr"                         
[200] "prior_b_MASS"                              
[201] "prior_sigma"                               
[202] "prior_sd_FISHID"                           
[203] "prior_cor_FISHID"                          
[204] "lprior"                                    
[205] "lp__"                                      
[206] "accept_stat__"                             
[207] "stepsize__"                                
[208] "treedepth__"                               
[209] "n_leapfrog__"                              
[210] "divergent__"                               
[211] "energy__"

pars <- norin.brm4 |> get_variables()
pars <- str_extract(pars, "^b_.*|^sigma$|^sd.*") |> na.omit()

norin.brm4$fit |>
  stan_trace(pars = pars)

The chains appear well mixed and very similar

stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

norin.brm4$fit |>
  stan_ac(pars = pars)

There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

norin.brm4$fit |> stan_rhat()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

norin.brm4$fit |> stan_ess()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.

norin.brm4$fit |>
  stan_dens(separate_chains = TRUE, pars = pars)

The ggmean package also as a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

norin.ggs <- norin.brm4 |> ggs(burnin = FALSE, inc_warmup = FALSE)

Warning in custom.sort(D$Parameter): NAs introduced by coercion

norin.ggs |> ggs_traceplot()

The chains appear well mixed and very similar

gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(norin.ggs)

There is no evidence of autocorrelation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(norin.ggs, scaling = 1.01)

All Rhat values are below 1.01, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

ggs_effective(norin.ggs)

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
dplyr 1.1.0.
ℹ Please use `reframe()` instead.
ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
  always returns an ungrouped data frame and adjust accordingly.
ℹ The deprecated feature was likely used in the ggmcmc package.
  Please report the issue at <https://github.com/xfim/ggmcmc/issues/>.

Ratios all very high.

More diagnostics

ggs_crosscorrelation(norin.ggs)

ggs_grb(norin.ggs)

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

See list of available diagnostics by name

available_ppc()

bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped

dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

norin.brm4 |> pp_check(type = "dens_overlay", ndraws = 100)

The model draws appear to be consistent with the observed data.

error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. Note, this is not really that useful for models that involve a binomial response

norin.brm4 |> pp_check(type = "error_scatter_avg")

Using all posterior draws for ppc type 'error_scatter_avg' by default.

This is not really interpretable

intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

norin.brm4 |> pp_check(group = "TRIAL", type = "intervals")

Using all posterior draws for ppc type 'intervals' by default.

Warning: The following arguments were unrecognized and ignored: group

norin.brm3 |> pp_check(group = "TRIAL", type = "intervals_grouped")

Using all posterior draws for ppc type 'intervals_grouped' by default.

norin.brm3 |> pp_check(group = "TRIAL", type = "violin_grouped")

Using all posterior draws for ppc type 'violin_grouped' by default.

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(norin.brm2)

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

simulated (predicted) responses associated with each observation.
observed values
fitted (predicted) responses (averaged) associated with each observation

preds <- norin.brm4 |> posterior_predict(ndraws = 250, summary = FALSE)
norin.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = norin$CHANGE,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = FALSE
)
plot(norin.resids, quantreg = FALSE)

norin.resids |> testDispersion()


    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.2998, p-value < 2.2e-16
alternative hypothesis: two.sided

norin.resids <- make_brms_dharma_res(norin.brm4, integerResponse = FALSE)
wrap_elements(~ testUniformity(norin.resids)) +
  wrap_elements(~ plotResiduals(norin.resids, form = factor(rep(1, nrow(norin))))) +
  wrap_elements(~ plotResiduals(norin.resids)) +
  wrap_elements(~ testDispersion(norin.resids))

Conclusions:

the simulated residuals do suggest some issues with the residuals
the model appears to be under-dispersed
there is evidence of a lack of fit.

7 Partial effects plots

brms

norin.brm4 |>
  conditional_effects() |>
  plot(ask = FALSE, points = TRUE, plot = FALSE) |>
  wrap_plots()

norin.brm4 |>
  ggpredict(terms = c("SMR_contr", "TRIAL")) |>
  plot(show_data = TRUE)

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

norin.brm4 |>
  ggemmeans(~ SMR_contr * TRIAL) |>
  plot(show_data = TRUE)

Warning in stats::qt(ci, df = dof): NaNs produced

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

It is not really possible to do this via the fitted draws as it would not be marginalising over MASS or FISHID.

8 Model investigation

brms

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

norin.brm4 |> summary()

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + MASS 
   Data: norin (Number of observations: 180) 
  Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
         total post-warmup draws = 2400

Multilevel Hyperparameters:
~FISHID (Number of levels: 60) 
                                   Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                         11.77      2.94     5.72    17.45 1.01
sd(TRIALHypoxia)                       6.31      3.72     0.39    14.59 1.00
sd(TRIALLowSalinity)                  28.09      3.90    20.80    36.43 1.00
cor(Intercept,TRIALHypoxia)           -0.33      0.38    -0.88     0.62 1.00
cor(Intercept,TRIALLowSalinity)        0.08      0.24    -0.34     0.61 1.00
cor(TRIALHypoxia,TRIALLowSalinity)    -0.48      0.37    -0.95     0.52 1.00
                                   Bulk_ESS Tail_ESS
sd(Intercept)                           729     1248
sd(TRIALHypoxia)                        725      672
sd(TRIALLowSalinity)                   1501     1742
cor(Intercept,TRIALHypoxia)            1857     2050
cor(Intercept,TRIALLowSalinity)        1368     1600
cor(TRIALHypoxia,TRIALLowSalinity)      611      519

Regression Coefficients:
                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                    147.10     22.28   103.34   190.74 1.00     2420
TRIALHypoxia                -107.54     21.45  -149.71   -64.56 1.00     2372
TRIALLowSalinity             -72.63     32.21  -134.70    -5.92 1.00     2262
SMR_contr                    -19.59      3.68   -26.60   -12.29 1.00     2433
MASS                           0.12      0.24    -0.34     0.59 1.00     2409
TRIALHypoxia:SMR_contr        10.96      4.14     2.86    18.97 1.00     2370
TRIALLowSalinity:SMR_contr     5.95      6.19    -6.63    17.75 1.00     2278
                           Tail_ESS
Intercept                      2292
TRIALHypoxia                   2227
TRIALLowSalinity               2002
SMR_contr                      1999
MASS                           2288
TRIALHypoxia:SMR_contr         2383
TRIALLowSalinity:SMR_contr     2159

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    13.37      1.94     9.50    17.21 1.01      549      517

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

the intercept indicates that the average metabolic change associated with the High Temperature trial (at the mean SMR_contr and MASS) across all fish was 147.1.
for the same SMR_contr and MASS, the average metabolic change was 107.54 less in the Hypoxia trial and 72.63 less in the Salinity trial than the High Temperature trial.
in the High Temperature trial, every 1 unit increase in SMR_contr was associated with a decline in metabolic change of round(-1*norin.sum$fixed[4,1], 2) units.
in the High Temperature trial (and mean SMR_contr), every 1 unit increase in MASS was associated with an increase in metabolic change of round(-1*norin.sum$fixed[5,1], 2) units.
the rate of change in metabolic change associated with the Hypoxia trial was round(-1*norin.sum$fixed[6,1], 2) units higher than that of the High Temperature trial and round(-1*norin.sum$fixed[7,1], 2) units higher in the Low Salinity trial.
the variance in intercepts across all fish is 11.77
the variance in responses between High Temperature and Hypoxia trials across all fish is 6.31.
the variance in responses between High Temperature and Low Salinity trials across all fish is 28.09.
the scale of variance between fish within a treatment
(sigma, 13.37) is substantially higher than that both within the High Temperature trial and between the High Temperature and Hypoxia trials, yet lower than that between High Temperature and Low Salinity trials.

norin.brm4 |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    rhat,
    ess_bulk,
    ess_tail
  )

# A tibble: 205 × 7
   variable                      median    lower   upper  rhat ess_bulk ess_tail
   <chr>                          <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 b_Intercept                  146.     1.04e+2 191.    1.00     2420.    2292.
 2 b_TRIALHypoxia              -108.    -1.48e+2 -64.1   1.00     2372.    2227.
 3 b_TRIALLowSalinity           -72.9   -1.37e+2 -10.2   1.000    2262.    2002.
 4 b_SMR_contr                  -19.7   -2.67e+1 -12.4   1.000    2433.    1999.
 5 b_MASS                         0.118 -3.51e-1   0.573 1.000    2408.    2288.
 6 b_TRIALHypoxia:SMR_contr      11.0    2.84e+0  18.9   1.00     2370.    2383.
 7 b_TRIALLowSalinity:SMR_con…    5.89  -6.37e+0  17.8   0.999    2279.    2159.
 8 sd_FISHID__Intercept          11.8    5.64e+0  17.3   1.01      729.    1248.
 9 sd_FISHID__TRIALHypoxia        6.07   2.69e-3  12.8   1.00      725.     672.
10 sd_FISHID__TRIALLowSalinity   28.0    2.03e+1  35.5   1.00     1500.    1742.
# ℹ 195 more rows

## or if you want to exclude some parameters
norin.brm4 |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    rhat,
    ess_bulk,
    ess_tail
  ) |>
  filter(str_detect(variable, "prior|^r_|^lp__", negate = TRUE))

# A tibble: 15 × 7
   variable                      median    lower   upper  rhat ess_bulk ess_tail
   <chr>                          <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 b_Intercept                  1.46e+2  1.04e+2 191.    1.00     2420.    2292.
 2 b_TRIALHypoxia              -1.08e+2 -1.48e+2 -64.1   1.00     2372.    2227.
 3 b_TRIALLowSalinity          -7.29e+1 -1.37e+2 -10.2   1.000    2262.    2002.
 4 b_SMR_contr                 -1.97e+1 -2.67e+1 -12.4   1.000    2433.    1999.
 5 b_MASS                       1.18e-1 -3.51e-1   0.573 1.000    2408.    2288.
 6 b_TRIALHypoxia:SMR_contr     1.10e+1  2.84e+0  18.9   1.00     2370.    2383.
 7 b_TRIALLowSalinity:SMR_con…  5.89e+0 -6.37e+0  17.8   0.999    2279.    2159.
 8 sd_FISHID__Intercept         1.18e+1  5.64e+0  17.3   1.01      729.    1248.
 9 sd_FISHID__TRIALHypoxia      6.07e+0  2.69e-3  12.8   1.00      725.     672.
10 sd_FISHID__TRIALLowSalinity  2.80e+1  2.03e+1  35.5   1.00     1500.    1742.
11 cor_FISHID__Intercept__TRI… -3.97e-1 -9.58e-1   0.446 1.00     1856.    2050.
12 cor_FISHID__Intercept__TRI…  4.94e-2 -3.61e-1   0.574 1.00     1369.    1600.
13 cor_FISHID__TRIALHypoxia__… -5.50e-1 -9.99e-1   0.249 1.00      611.     519.
14 sigma                        1.34e+1  9.83e+0  17.4   1.01      548.     517.
15 Intercept                    1.97e+1  1.54e+1  23.5   0.999    2504.    2228.

norin.brm4 |> as_draws_df()

# A draws_df: 800 iterations, 3 chains, and 205 variables
   b_Intercept b_TRIALHypoxia b_TRIALLowSalinity b_SMR_contr  b_MASS
1          139           -153                -73         -20  0.3392
2          132           -103                -73         -17  0.2897
3          135           -107                -39         -19  0.3357
4          169           -140                -81         -23 -0.0097
5          144           -127                 -3         -19  0.1346
6          154           -104                  1         -22  0.2451
7          178           -105                -48         -24 -0.0865
8          155           -105                -56         -21  0.0615
9          132            -80                -62         -15  0.0266
10         149           -106               -109         -21  0.2686
   b_TRIALHypoxia:SMR_contr b_TRIALLowSalinity:SMR_contr sd_FISHID__Intercept
1                      19.9                          7.1                 10.8
2                       9.8                          6.5                 10.8
3                      10.0                         -1.2                 10.0
4                      17.0                          7.6                 15.0
5                      14.9                         -7.3                 12.0
6                      10.3                         -9.0                 11.5
7                      10.3                          1.8                 10.3
8                      10.7                          3.1                 13.0
9                       4.8                          4.2                 15.1
10                     10.3                         12.5                  8.3
# ... with 2390 more draws, and 197 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}

norin.brm4 |>
  as_draws_df() |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    rhat,
    ess_bulk,
    ess_tail
  )

# A tibble: 205 × 7
   variable                      median    lower   upper  rhat ess_bulk ess_tail
   <chr>                          <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 b_Intercept                  146.     1.04e+2 191.    1.00     2420.    2292.
 2 b_TRIALHypoxia              -108.    -1.48e+2 -64.1   1.00     2372.    2227.
 3 b_TRIALLowSalinity           -72.9   -1.37e+2 -10.2   1.000    2262.    2002.
 4 b_SMR_contr                  -19.7   -2.67e+1 -12.4   1.000    2433.    1999.
 5 b_MASS                         0.118 -3.51e-1   0.573 1.000    2408.    2288.
 6 b_TRIALHypoxia:SMR_contr      11.0    2.84e+0  18.9   1.00     2370.    2383.
 7 b_TRIALLowSalinity:SMR_con…    5.89  -6.37e+0  17.8   0.999    2279.    2159.
 8 sd_FISHID__Intercept          11.8    5.64e+0  17.3   1.01      729.    1248.
 9 sd_FISHID__TRIALHypoxia        6.07   2.69e-3  12.8   1.00      725.     672.
10 sd_FISHID__TRIALLowSalinity   28.0    2.03e+1  35.5   1.00     1500.    1742.
# ℹ 195 more rows

## or if you want to exclude some parameters
norin.brm4 |>
  as_draws_df() |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    rhat,
    ess_bulk,
    ess_tail
  ) |>
  filter(str_detect(variable, "prior|^r_|^lp__", negate = TRUE))

# A tibble: 15 × 7
   variable                      median    lower   upper  rhat ess_bulk ess_tail
   <chr>                          <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 b_Intercept                  1.46e+2  1.04e+2 191.    1.00     2420.    2292.
 2 b_TRIALHypoxia              -1.08e+2 -1.48e+2 -64.1   1.00     2372.    2227.
 3 b_TRIALLowSalinity          -7.29e+1 -1.37e+2 -10.2   1.000    2262.    2002.
 4 b_SMR_contr                 -1.97e+1 -2.67e+1 -12.4   1.000    2433.    1999.
 5 b_MASS                       1.18e-1 -3.51e-1   0.573 1.000    2408.    2288.
 6 b_TRIALHypoxia:SMR_contr     1.10e+1  2.84e+0  18.9   1.00     2370.    2383.
 7 b_TRIALLowSalinity:SMR_con…  5.89e+0 -6.37e+0  17.8   0.999    2279.    2159.
 8 sd_FISHID__Intercept         1.18e+1  5.64e+0  17.3   1.01      729.    1248.
 9 sd_FISHID__TRIALHypoxia      6.07e+0  2.69e-3  12.8   1.00      725.     672.
10 sd_FISHID__TRIALLowSalinity  2.80e+1  2.03e+1  35.5   1.00     1500.    1742.
11 cor_FISHID__Intercept__TRI… -3.97e-1 -9.58e-1   0.446 1.00     1856.    2050.
12 cor_FISHID__Intercept__TRI…  4.94e-2 -3.61e-1   0.574 1.00     1369.    1600.
13 cor_FISHID__TRIALHypoxia__… -5.50e-1 -9.99e-1   0.249 1.00      611.     519.
14 sigma                        1.34e+1  9.83e+0  17.4   1.01      548.     517.
15 Intercept                    1.97e+1  1.54e+1  23.5   0.999    2504.    2228.

norin.brm4$fit |>
  tidyMCMC(
    estimate.method = "median",
    conf.int = TRUE, conf.method = "HPDinterval",
    rhat = TRUE, ess = TRUE
  )

# A tibble: 204 × 7
   term                        estimate std.error conf.low conf.high  rhat   ess
   <chr>                          <dbl>     <dbl>    <dbl>     <dbl> <dbl> <int>
 1 b_Intercept                  147.       22.3    1.04e+2   191.    0.999  2412
 2 b_TRIALHypoxia              -108.       21.4   -1.48e+2   -64.1   1.000  2366
 3 b_TRIALLowSalinity           -72.6      32.2   -1.37e+2   -10.2   0.999  2234
 4 b_SMR_contr                  -19.6       3.68  -2.67e+1   -12.4   0.999  2427
 5 b_MASS                         0.121     0.240 -3.51e-1     0.573 1.000  2414
 6 b_TRIALHypoxia:SMR_contr      11.0       4.14   2.84e+0    18.9   1.000  2358
 7 b_TRIALLowSalinity:SMR_con…    5.95      6.19  -6.37e+0    17.8   0.999  2245
 8 sd_FISHID__Intercept          11.8       2.94   5.64e+0    17.3   1.01    633
 9 sd_FISHID__TRIALHypoxia        6.31      3.72   2.69e-3    12.8   1.00    704
10 sd_FISHID__TRIALLowSalinity   28.1       3.90   2.03e+1    35.5   1.00   1499
# ℹ 194 more rows

Conclusions:

see above

norin.brm4 |> get_variables()

  [1] "b_Intercept"                               
  [2] "b_TRIALHypoxia"                            
  [3] "b_TRIALLowSalinity"                        
  [4] "b_SMR_contr"                               
  [5] "b_MASS"                                    
  [6] "b_TRIALHypoxia:SMR_contr"                  
  [7] "b_TRIALLowSalinity:SMR_contr"              
  [8] "sd_FISHID__Intercept"                      
  [9] "sd_FISHID__TRIALHypoxia"                   
 [10] "sd_FISHID__TRIALLowSalinity"               
 [11] "cor_FISHID__Intercept__TRIALHypoxia"       
 [12] "cor_FISHID__Intercept__TRIALLowSalinity"   
 [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
 [14] "sigma"                                     
 [15] "Intercept"                                 
 [16] "r_FISHID[1,Intercept]"                     
 [17] "r_FISHID[2,Intercept]"                     
 [18] "r_FISHID[3,Intercept]"                     
 [19] "r_FISHID[4,Intercept]"                     
 [20] "r_FISHID[5,Intercept]"                     
 [21] "r_FISHID[6,Intercept]"                     
 [22] "r_FISHID[7,Intercept]"                     
 [23] "r_FISHID[8,Intercept]"                     
 [24] "r_FISHID[9,Intercept]"                     
 [25] "r_FISHID[10,Intercept]"                    
 [26] "r_FISHID[11,Intercept]"                    
 [27] "r_FISHID[12,Intercept]"                    
 [28] "r_FISHID[13,Intercept]"                    
 [29] "r_FISHID[14,Intercept]"                    
 [30] "r_FISHID[15,Intercept]"                    
 [31] "r_FISHID[16,Intercept]"                    
 [32] "r_FISHID[17,Intercept]"                    
 [33] "r_FISHID[18,Intercept]"                    
 [34] "r_FISHID[19,Intercept]"                    
 [35] "r_FISHID[20,Intercept]"                    
 [36] "r_FISHID[21,Intercept]"                    
 [37] "r_FISHID[22,Intercept]"                    
 [38] "r_FISHID[23,Intercept]"                    
 [39] "r_FISHID[24,Intercept]"                    
 [40] "r_FISHID[25,Intercept]"                    
 [41] "r_FISHID[26,Intercept]"                    
 [42] "r_FISHID[27,Intercept]"                    
 [43] "r_FISHID[28,Intercept]"                    
 [44] "r_FISHID[29,Intercept]"                    
 [45] "r_FISHID[30,Intercept]"                    
 [46] "r_FISHID[31,Intercept]"                    
 [47] "r_FISHID[32,Intercept]"                    
 [48] "r_FISHID[33,Intercept]"                    
 [49] "r_FISHID[34,Intercept]"                    
 [50] "r_FISHID[35,Intercept]"                    
 [51] "r_FISHID[36,Intercept]"                    
 [52] "r_FISHID[37,Intercept]"                    
 [53] "r_FISHID[38,Intercept]"                    
 [54] "r_FISHID[39,Intercept]"                    
 [55] "r_FISHID[40,Intercept]"                    
 [56] "r_FISHID[41,Intercept]"                    
 [57] "r_FISHID[42,Intercept]"                    
 [58] "r_FISHID[43,Intercept]"                    
 [59] "r_FISHID[44,Intercept]"                    
 [60] "r_FISHID[45,Intercept]"                    
 [61] "r_FISHID[46,Intercept]"                    
 [62] "r_FISHID[47,Intercept]"                    
 [63] "r_FISHID[48,Intercept]"                    
 [64] "r_FISHID[49,Intercept]"                    
 [65] "r_FISHID[50,Intercept]"                    
 [66] "r_FISHID[51,Intercept]"                    
 [67] "r_FISHID[52,Intercept]"                    
 [68] "r_FISHID[53,Intercept]"                    
 [69] "r_FISHID[54,Intercept]"                    
 [70] "r_FISHID[55,Intercept]"                    
 [71] "r_FISHID[56,Intercept]"                    
 [72] "r_FISHID[57,Intercept]"                    
 [73] "r_FISHID[58,Intercept]"                    
 [74] "r_FISHID[59,Intercept]"                    
 [75] "r_FISHID[60,Intercept]"                    
 [76] "r_FISHID[1,TRIALHypoxia]"                  
 [77] "r_FISHID[2,TRIALHypoxia]"                  
 [78] "r_FISHID[3,TRIALHypoxia]"                  
 [79] "r_FISHID[4,TRIALHypoxia]"                  
 [80] "r_FISHID[5,TRIALHypoxia]"                  
 [81] "r_FISHID[6,TRIALHypoxia]"                  
 [82] "r_FISHID[7,TRIALHypoxia]"                  
 [83] "r_FISHID[8,TRIALHypoxia]"                  
 [84] "r_FISHID[9,TRIALHypoxia]"                  
 [85] "r_FISHID[10,TRIALHypoxia]"                 
 [86] "r_FISHID[11,TRIALHypoxia]"                 
 [87] "r_FISHID[12,TRIALHypoxia]"                 
 [88] "r_FISHID[13,TRIALHypoxia]"                 
 [89] "r_FISHID[14,TRIALHypoxia]"                 
 [90] "r_FISHID[15,TRIALHypoxia]"                 
 [91] "r_FISHID[16,TRIALHypoxia]"                 
 [92] "r_FISHID[17,TRIALHypoxia]"                 
 [93] "r_FISHID[18,TRIALHypoxia]"                 
 [94] "r_FISHID[19,TRIALHypoxia]"                 
 [95] "r_FISHID[20,TRIALHypoxia]"                 
 [96] "r_FISHID[21,TRIALHypoxia]"                 
 [97] "r_FISHID[22,TRIALHypoxia]"                 
 [98] "r_FISHID[23,TRIALHypoxia]"                 
 [99] "r_FISHID[24,TRIALHypoxia]"                 
[100] "r_FISHID[25,TRIALHypoxia]"                 
[101] "r_FISHID[26,TRIALHypoxia]"                 
[102] "r_FISHID[27,TRIALHypoxia]"                 
[103] "r_FISHID[28,TRIALHypoxia]"                 
[104] "r_FISHID[29,TRIALHypoxia]"                 
[105] "r_FISHID[30,TRIALHypoxia]"                 
[106] "r_FISHID[31,TRIALHypoxia]"                 
[107] "r_FISHID[32,TRIALHypoxia]"                 
[108] "r_FISHID[33,TRIALHypoxia]"                 
[109] "r_FISHID[34,TRIALHypoxia]"                 
[110] "r_FISHID[35,TRIALHypoxia]"                 
[111] "r_FISHID[36,TRIALHypoxia]"                 
[112] "r_FISHID[37,TRIALHypoxia]"                 
[113] "r_FISHID[38,TRIALHypoxia]"                 
[114] "r_FISHID[39,TRIALHypoxia]"                 
[115] "r_FISHID[40,TRIALHypoxia]"                 
[116] "r_FISHID[41,TRIALHypoxia]"                 
[117] "r_FISHID[42,TRIALHypoxia]"                 
[118] "r_FISHID[43,TRIALHypoxia]"                 
[119] "r_FISHID[44,TRIALHypoxia]"                 
[120] "r_FISHID[45,TRIALHypoxia]"                 
[121] "r_FISHID[46,TRIALHypoxia]"                 
[122] "r_FISHID[47,TRIALHypoxia]"                 
[123] "r_FISHID[48,TRIALHypoxia]"                 
[124] "r_FISHID[49,TRIALHypoxia]"                 
[125] "r_FISHID[50,TRIALHypoxia]"                 
[126] "r_FISHID[51,TRIALHypoxia]"                 
[127] "r_FISHID[52,TRIALHypoxia]"                 
[128] "r_FISHID[53,TRIALHypoxia]"                 
[129] "r_FISHID[54,TRIALHypoxia]"                 
[130] "r_FISHID[55,TRIALHypoxia]"                 
[131] "r_FISHID[56,TRIALHypoxia]"                 
[132] "r_FISHID[57,TRIALHypoxia]"                 
[133] "r_FISHID[58,TRIALHypoxia]"                 
[134] "r_FISHID[59,TRIALHypoxia]"                 
[135] "r_FISHID[60,TRIALHypoxia]"                 
[136] "r_FISHID[1,TRIALLowSalinity]"              
[137] "r_FISHID[2,TRIALLowSalinity]"              
[138] "r_FISHID[3,TRIALLowSalinity]"              
[139] "r_FISHID[4,TRIALLowSalinity]"              
[140] "r_FISHID[5,TRIALLowSalinity]"              
[141] "r_FISHID[6,TRIALLowSalinity]"              
[142] "r_FISHID[7,TRIALLowSalinity]"              
[143] "r_FISHID[8,TRIALLowSalinity]"              
[144] "r_FISHID[9,TRIALLowSalinity]"              
[145] "r_FISHID[10,TRIALLowSalinity]"             
[146] "r_FISHID[11,TRIALLowSalinity]"             
[147] "r_FISHID[12,TRIALLowSalinity]"             
[148] "r_FISHID[13,TRIALLowSalinity]"             
[149] "r_FISHID[14,TRIALLowSalinity]"             
[150] "r_FISHID[15,TRIALLowSalinity]"             
[151] "r_FISHID[16,TRIALLowSalinity]"             
[152] "r_FISHID[17,TRIALLowSalinity]"             
[153] "r_FISHID[18,TRIALLowSalinity]"             
[154] "r_FISHID[19,TRIALLowSalinity]"             
[155] "r_FISHID[20,TRIALLowSalinity]"             
[156] "r_FISHID[21,TRIALLowSalinity]"             
[157] "r_FISHID[22,TRIALLowSalinity]"             
[158] "r_FISHID[23,TRIALLowSalinity]"             
[159] "r_FISHID[24,TRIALLowSalinity]"             
[160] "r_FISHID[25,TRIALLowSalinity]"             
[161] "r_FISHID[26,TRIALLowSalinity]"             
[162] "r_FISHID[27,TRIALLowSalinity]"             
[163] "r_FISHID[28,TRIALLowSalinity]"             
[164] "r_FISHID[29,TRIALLowSalinity]"             
[165] "r_FISHID[30,TRIALLowSalinity]"             
[166] "r_FISHID[31,TRIALLowSalinity]"             
[167] "r_FISHID[32,TRIALLowSalinity]"             
[168] "r_FISHID[33,TRIALLowSalinity]"             
[169] "r_FISHID[34,TRIALLowSalinity]"             
[170] "r_FISHID[35,TRIALLowSalinity]"             
[171] "r_FISHID[36,TRIALLowSalinity]"             
[172] "r_FISHID[37,TRIALLowSalinity]"             
[173] "r_FISHID[38,TRIALLowSalinity]"             
[174] "r_FISHID[39,TRIALLowSalinity]"             
[175] "r_FISHID[40,TRIALLowSalinity]"             
[176] "r_FISHID[41,TRIALLowSalinity]"             
[177] "r_FISHID[42,TRIALLowSalinity]"             
[178] "r_FISHID[43,TRIALLowSalinity]"             
[179] "r_FISHID[44,TRIALLowSalinity]"             
[180] "r_FISHID[45,TRIALLowSalinity]"             
[181] "r_FISHID[46,TRIALLowSalinity]"             
[182] "r_FISHID[47,TRIALLowSalinity]"             
[183] "r_FISHID[48,TRIALLowSalinity]"             
[184] "r_FISHID[49,TRIALLowSalinity]"             
[185] "r_FISHID[50,TRIALLowSalinity]"             
[186] "r_FISHID[51,TRIALLowSalinity]"             
[187] "r_FISHID[52,TRIALLowSalinity]"             
[188] "r_FISHID[53,TRIALLowSalinity]"             
[189] "r_FISHID[54,TRIALLowSalinity]"             
[190] "r_FISHID[55,TRIALLowSalinity]"             
[191] "r_FISHID[56,TRIALLowSalinity]"             
[192] "r_FISHID[57,TRIALLowSalinity]"             
[193] "r_FISHID[58,TRIALLowSalinity]"             
[194] "r_FISHID[59,TRIALLowSalinity]"             
[195] "r_FISHID[60,TRIALLowSalinity]"             
[196] "prior_Intercept"                           
[197] "prior_b_TRIALHypoxia"                      
[198] "prior_b_TRIALLowSalinity"                  
[199] "prior_b_SMR_contr"                         
[200] "prior_b_MASS"                              
[201] "prior_sigma"                               
[202] "prior_sd_FISHID"                           
[203] "prior_cor_FISHID"                          
[204] "lprior"                                    
[205] "lp__"                                      
[206] "accept_stat__"                             
[207] "stepsize__"                                
[208] "treedepth__"                               
[209] "n_leapfrog__"                              
[210] "divergent__"                               
[211] "energy__"

norin.draw <- norin.brm4 |>
  gather_draws(`^b.Intercept$|b_.*|sd_.*|sigma`, regex = TRUE)
norin.draw

# A tibble: 26,400 × 5
# Groups:   .variable [11]
   .chain .iteration .draw .variable   .value
    <int>      <int> <int> <chr>        <dbl>
 1      1          1     1 b_Intercept   139.
 2      1          2     2 b_Intercept   132.
 3      1          3     3 b_Intercept   135.
 4      1          4     4 b_Intercept   169.
 5      1          5     5 b_Intercept   144.
 6      1          6     6 b_Intercept   154.
 7      1          7     7 b_Intercept   178.
 8      1          8     8 b_Intercept   155.
 9      1          9     9 b_Intercept   132.
10      1         10    10 b_Intercept   149.
# ℹ 26,390 more rows

We can then summarise this

norin.draw |> median_hdci()

# A tibble: 11 × 7
   .variable                     .value   .lower  .upper .width .point .interval
   <chr>                          <dbl>    <dbl>   <dbl>  <dbl> <chr>  <chr>    
 1 b_Intercept                  146.     1.01e+2 189.      0.95 median hdci     
 2 b_MASS                         0.118 -3.51e-1   0.573   0.95 median hdci     
 3 b_SMR_contr                  -19.7   -2.67e+1 -12.4     0.95 median hdci     
 4 b_TRIALHypoxia              -108.    -1.53e+2 -68.3     0.95 median hdci     
 5 b_TRIALHypoxia:SMR_contr      11.0    2.51e+0  18.6     0.95 median hdci     
 6 b_TRIALLowSalinity           -72.9   -1.37e+2 -10.2     0.95 median hdci     
 7 b_TRIALLowSalinity:SMR_con…    5.89  -6.65e+0  17.7     0.95 median hdci     
 8 sd_FISHID__Intercept          11.8    6.08e+0  17.8     0.95 median hdci     
 9 sd_FISHID__TRIALHypoxia        6.07   2.69e-3  12.8     0.95 median hdci     
10 sd_FISHID__TRIALLowSalinity   28.0    2.03e+1  35.5     0.95 median hdci     
11 sigma                         13.4    9.83e+0  17.4     0.95 median hdci

norin.brm4 |>
  gather_draws(`b_Intercept|b_.*`, regex = TRUE) |>
  ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  stat_slab(aes(
    x = .value, y = .variable,
    fill = stat(ggdist::cut_cdf_qi(cdf,
      .width = c(0.5, 0.8, 0.95),
      labels = scales::percent_format()
    ))
  ), color = "black") +
  scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)

Warning: `stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95), labels =
scales::percent_format()))` was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95),
  labels = scales::percent_format()))` instead.

norin.brm4 |>
  gather_draws(`b_Intercept|b_.*`, regex = TRUE) |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  facet_wrap(~.variable, scales = "free")

norin.brm4 |>
  gather_draws(`b_Intercept|b_.*`, regex = TRUE) |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  theme_classic()

norin.brm4$fit |> plot(type = "intervals")

'pars' not specified. Showing first 10 parameters by default.

ci_level: 0.8 (80% intervals)

outer_level: 0.95 (95% intervals)

norin.brm4 |>
  gather_draws(`^b_.*`, regex = TRUE) |>
  filter(.variable != "b_Intercept") |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  facet_wrap(~.variable, scales = "free")

norin.brm4 |>
  gather_draws(`^b_.*`, regex = TRUE) |>
  filter(.variable != "b_Intercept") |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  geom_vline(xintercept = 0, linetype = "dashed")

norin.brm4 |>
  gather_draws(`^b_.*`, regex = TRUE) |>
  filter(.variable != "b_Intercept") |>
  ggplot() +
  geom_density_ridges(aes(x = .value, y = .variable), alpha = 0.4) +
  geom_vline(xintercept = 0, linetype = "dashed")

Picking joint bandwidth of 2.12

## Or in colour
norin.brm4 |>
  gather_draws(`^b_.*`, regex = TRUE) |>
  filter(.variable != "b_Intercept") |>
  ggplot() +
  geom_density_ridges_gradient(
    aes(
      x = exp(.value),
      y = .variable,
      fill = stat(x)
    ),
    alpha = 0.4, colour = "white",
    quantile_lines = TRUE,
    quantiles = c(0.025, 0.975)
  ) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  scale_x_continuous(trans = scales::log2_trans()) +
  scale_fill_viridis_c(option = "C")

Picking joint bandwidth of 3.06

norin.brm4 |> tidy_draws()

# A tibble: 2,400 × 214
   .chain .iteration .draw b_Intercept b_TRIALHypoxia b_TRIALLowSalinity
    <int>      <int> <int>       <dbl>          <dbl>              <dbl>
 1      1          1     1        139.         -153.              -73.1 
 2      1          2     2        132.         -103.              -73.4 
 3      1          3     3        135.         -107.              -39.0 
 4      1          4     4        169.         -140.              -81.1 
 5      1          5     5        144.         -127.               -3.04
 6      1          6     6        154.         -104.                1.02
 7      1          7     7        178.         -105.              -47.8 
 8      1          8     8        155.         -105.              -55.8 
 9      1          9     9        132.          -80.0             -61.9 
10      1         10    10        149.         -106.             -109.  
# ℹ 2,390 more rows
# ℹ 208 more variables: b_SMR_contr <dbl>, b_MASS <dbl>,
#   `b_TRIALHypoxia:SMR_contr` <dbl>, `b_TRIALLowSalinity:SMR_contr` <dbl>,
#   sd_FISHID__Intercept <dbl>, sd_FISHID__TRIALHypoxia <dbl>,
#   sd_FISHID__TRIALLowSalinity <dbl>,
#   cor_FISHID__Intercept__TRIALHypoxia <dbl>,
#   cor_FISHID__Intercept__TRIALLowSalinity <dbl>, …

norin.brm4 |> spread_draws(`.*Intercept.*|^b_.*`, regex = TRUE)

# A tibble: 2,400 × 75
   .chain .iteration .draw b_Intercept b_TRIALHypoxia b_TRIALLowSalinity
    <int>      <int> <int>       <dbl>          <dbl>              <dbl>
 1      1          1     1        139.         -153.              -73.1 
 2      1          2     2        132.         -103.              -73.4 
 3      1          3     3        135.         -107.              -39.0 
 4      1          4     4        169.         -140.              -81.1 
 5      1          5     5        144.         -127.               -3.04
 6      1          6     6        154.         -104.                1.02
 7      1          7     7        178.         -105.              -47.8 
 8      1          8     8        155.         -105.              -55.8 
 9      1          9     9        132.          -80.0             -61.9 
10      1         10    10        149.         -106.             -109.  
# ℹ 2,390 more rows
# ℹ 69 more variables: b_SMR_contr <dbl>, b_MASS <dbl>,
#   `b_TRIALHypoxia:SMR_contr` <dbl>, `b_TRIALLowSalinity:SMR_contr` <dbl>,
#   sd_FISHID__Intercept <dbl>, cor_FISHID__Intercept__TRIALHypoxia <dbl>,
#   cor_FISHID__Intercept__TRIALLowSalinity <dbl>, Intercept <dbl>,
#   `r_FISHID[1,Intercept]` <dbl>, `r_FISHID[2,Intercept]` <dbl>,
#   `r_FISHID[3,Intercept]` <dbl>, `r_FISHID[4,Intercept]` <dbl>, …

norin.brm4 |>
  posterior_samples() |>
  as_tibble()

Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.

# A tibble: 2,400 × 205
   b_Intercept b_TRIALHypoxia b_TRIALLowSalinity b_SMR_contr   b_MASS
         <dbl>          <dbl>              <dbl>       <dbl>    <dbl>
 1        139.         -153.              -73.1        -19.9  0.339  
 2        132.         -103.              -73.4        -17.1  0.290  
 3        135.         -107.              -39.0        -18.5  0.336  
 4        169.         -140.              -81.1        -22.7 -0.00966
 5        144.         -127.               -3.04       -19.4  0.135  
 6        154.         -104.                1.02       -22.4  0.245  
 7        178.         -105.              -47.8        -23.9 -0.0865 
 8        155.         -105.              -55.8        -20.6  0.0615 
 9        132.          -80.0             -61.9        -15.1  0.0266 
10        149.         -106.             -109.         -21.3  0.269  
# ℹ 2,390 more rows
# ℹ 200 more variables: `b_TRIALHypoxia:SMR_contr` <dbl>,
#   `b_TRIALLowSalinity:SMR_contr` <dbl>, sd_FISHID__Intercept <dbl>,
#   sd_FISHID__TRIALHypoxia <dbl>, sd_FISHID__TRIALLowSalinity <dbl>,
#   cor_FISHID__Intercept__TRIALHypoxia <dbl>,
#   cor_FISHID__Intercept__TRIALLowSalinity <dbl>,
#   cor_FISHID__TRIALHypoxia__TRIALLowSalinity <dbl>, sigma <dbl>, …

norin.brm4 |>
  bayes_R2(re.form = NA, summary = FALSE) |>
  median_hdci()

          y      ymin      ymax .width .point .interval
1 0.4980791 0.4431186 0.5484889   0.95 median      hdci

norin.brm4 |>
  bayes_R2(re.form = ~ (1 | FISHID), summary = FALSE) |>
  median_hdci()

          y      ymin      ymax .width .point .interval
1 0.6128526 0.5304293 0.6687602   0.95 median      hdci

norin.brm4 |>
  bayes_R2(re.form = ~ (TRIAL | FISHID), summary = FALSE) |>
  median_hdci()

          y      ymin      ymax .width .point .interval
1 0.8459696 0.7474859 0.9157748   0.95 median      hdci

Region of Practical Equivalence

0.1 * sd(norin$CHANGE)

[1] 3.38454

norin.brm4 |> rope(range = c(-3.38, 3.38))

# Proportion of samples inside the ROPE [-3.38, 3.38]:

Parameter                  | inside ROPE
----------------------------------------
Intercept                  |      0.00 %
TRIALHypoxia               |      0.00 %
TRIALLowSalinity           |      0.00 %
SMR_contr                  |      0.00 %
MASS                       |    100.00 %
TRIALHypoxia:SMR_contr     |      0.70 %
TRIALLowSalinity:SMR_contr |     27.76 %

rope(norin.brm4, range = c(-3.38, 3.38)) |> plot()

norin.brm4 |> modelsummary(
  statistic = c("conf.low", "conf.high"),
  shape = term ~ statistic,
  exponentiate = FALSE
)

Warning: 
`modelsummary` uses the `performance` package to extract goodness-of-fit
statistics from models of this class. You can specify the statistics you wish
to compute by supplying a `metrics` argument to `modelsummary`, which will then
push it forward to `performance`. Acceptable values are: "all", "common",
"none", or a character vector of metrics names. For example: `modelsummary(mod,
metrics = c("RMSE", "R2")` Note that some metrics are computationally
expensive. See `?performance::performance` for details.
 This warning appears once per session.

	(1)
	Est.	2.5 %	97.5 %
b_Intercept	146.371	103.337	190.743
b_TRIALHypoxia	-107.503	-149.708	-64.556
b_TRIALLowSalinity	-72.901	-134.697	-5.921
b_SMR_contr	-19.665	-26.600	-12.287
b_MASS	0.118	-0.340	0.590
b_TRIALHypoxia × SMR_contr	10.981	2.857	18.973
b_TRIALLowSalinity × SMR_contr	5.895	-6.628	17.750
sd_FISHID__Intercept	11.770	5.723	17.449
sd_FISHID__TRIALHypoxia	6.071	0.392	14.586
sd_FISHID__TRIALLowSalinity	27.978	20.799	36.427
cor_FISHID__Intercept__TRIALHypoxia	-0.397	-0.881	0.616
cor_FISHID__Intercept__TRIALLowSalinity	0.049	-0.336	0.613
cor_FISHID__TRIALHypoxia__TRIALLowSalinity	-0.550	-0.952	0.519
sigma	13.382	9.501	17.208
Num.Obs.	180
R2	0.846
R2 Adj.	0.750
R2 Marg.	0.498
ELPD	-780.1
ELPD s.e.	9.2
LOOIC	1560.3
LOOIC s.e.	18.4
WAIC	1526.5
RMSE	9.26
r2.adjusted.marginal	0.308806833536584

norin.brm4 |> modelplot(exponentiate = FALSE)

9 Predictions / further analyses

brms

norin.brm4 |>
  emtrends(~TRIAL, var = "SMR_contr") |>
  pairs()

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

 contrast                      estimate lower.HPD upper.HPD
 HighTemperature - Hypoxia       -10.98    -18.94     -2.84
 HighTemperature - LowSalinity    -5.89    -17.83      6.37
 Hypoxia - LowSalinity             4.85     -8.74     17.46

Point estimate displayed: median 
HPD interval probability: 0.95

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

Conclusions:

the partial slope associated with control SMR during High Temperature was found to be 10.98 units more negative than during Hypoxia (although not significant). Note, we already had this partial slope in the original summary table.
the partial slope associated with control SMR during High Temperature was found to be 5.89 units more negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
the partial slope associated with control SMR during Hypoxia was found to be 4.85 units less negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.

We will estimated the pairwise differences between each of the three trials at the minimum, mean and maximum control SMR.

norin.grid <- with(norin, list(SMR_contr = Hmisc::smean.sdl(SMR_contr)))
norin.grid

$SMR_contr
    Mean    Lower    Upper 
5.178857 3.926671 6.431044

norin.brm4 |>
  emmeans(~ TRIAL | SMR_contr, at = norin.grid) |>
  pairs()

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

SMR_contr = 3.93:
 contrast                      estimate lower.HPD upper.HPD
 HighTemperature - Hypoxia        64.32      53.7     76.28
 HighTemperature - LowSalinity    49.41      33.1     66.86
 Hypoxia - LowSalinity           -15.09     -33.4      4.29

SMR_contr = 5.18:
 contrast                      estimate lower.HPD upper.HPD
 HighTemperature - Hypoxia        50.90      45.0     56.15
 HighTemperature - LowSalinity    41.71      33.3     50.35
 Hypoxia - LowSalinity            -8.89     -18.4      1.43

SMR_contr = 6.43:
 contrast                      estimate lower.HPD upper.HPD
 HighTemperature - Hypoxia        36.97      25.8     49.29
 HighTemperature - LowSalinity    34.31      16.2     51.17
 Hypoxia - LowSalinity            -2.86     -23.2     15.96

Point estimate displayed: median 
HPD interval probability: 0.95

norin.brm4 |>
  emmeans(~ TRIAL | SMR_contr, at = norin.grid) |>
  pairs() |>
  gather_emmeans_draws() |>
  median_hdci()

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

# A tibble: 9 × 8
  contrast                SMR_contr .value .lower .upper .width .point .interval
  <fct>                       <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 HighTemperature - Hypo…      3.93  64.3    53.7  76.3    0.95 median hdci     
2 HighTemperature - Hypo…      5.18  50.9    45.0  56.1    0.95 median hdci     
3 HighTemperature - Hypo…      6.43  37.0    25.8  49.3    0.95 median hdci     
4 HighTemperature - LowS…      3.93  49.4    33.1  66.9    0.95 median hdci     
5 HighTemperature - LowS…      5.18  41.7    33.3  50.4    0.95 median hdci     
6 HighTemperature - LowS…      6.43  34.3    16.5  51.7    0.95 median hdci     
7 Hypoxia - LowSalinity        3.93 -15.1   -33.6   4.15   0.95 median hdci     
8 Hypoxia - LowSalinity        5.18  -8.89  -18.6   1.33   0.95 median hdci     
9 Hypoxia - LowSalinity        6.43  -2.86  -23.2  16.0    0.95 median hdci

norin.brm4 |>
  emmeans(~ TRIAL | SMR_contr, at = norin.grid) |>
  pairs() |>
  tidy_draws() |>
  summarise_draws(median)

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

# A tibble: 9 × 2
  variable                                                          median
  <chr>                                                              <dbl>
1 contrast HighTemperature - Hypoxia SMR_contr 5.17885746916438      50.9 
2 contrast HighTemperature - LowSalinity SMR_contr 5.17885746916438  41.7 
3 contrast Hypoxia - LowSalinity SMR_contr 5.17885746916438          -8.89
4 contrast HighTemperature - Hypoxia SMR_contr 3.92667075536741      64.3 
5 contrast HighTemperature - LowSalinity SMR_contr 3.92667075536741  49.4 
6 contrast Hypoxia - LowSalinity SMR_contr 3.92667075536741         -15.1 
7 contrast HighTemperature - Hypoxia SMR_contr 6.43104418296134      37.0 
8 contrast HighTemperature - LowSalinity SMR_contr 6.43104418296134  34.3 
9 contrast Hypoxia - LowSalinity SMR_contr 6.43104418296134          -2.86

Conclusions:

as control SMR increases, the magnitude of differences between the trials becomes reduced.
the High Temperature trial resulted in the greatest positive change in SMR.
the change in SMR was not found to be different between the Hypoxia and Low Salinity trials.

Of course, as with other analyses, it is possible to extract the entire posterior and use it to calculate exceedence probabilities etc.

norin.em <- norin.brm4 |>
  emmeans(~ TRIAL | SMR_contr, at = norin.grid) |>
  pairs() |>
  gather_emmeans_draws() |>
  mutate(Fit = .value)

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

norin.em

# A tibble: 21,600 × 7
# Groups:   contrast, SMR_contr [9]
   contrast                  SMR_contr .chain .iteration .draw .value   Fit
   <fct>                         <dbl>  <int>      <int> <int>  <dbl> <dbl>
 1 HighTemperature - Hypoxia      5.18     NA         NA     1   49.6  49.6
 2 HighTemperature - Hypoxia      5.18     NA         NA     2   52.5  52.5
 3 HighTemperature - Hypoxia      5.18     NA         NA     3   55.3  55.3
 4 HighTemperature - Hypoxia      5.18     NA         NA     4   52.4  52.4
 5 HighTemperature - Hypoxia      5.18     NA         NA     5   49.7  49.7
 6 HighTemperature - Hypoxia      5.18     NA         NA     6   50.4  50.4
 7 HighTemperature - Hypoxia      5.18     NA         NA     7   51.3  51.3
 8 HighTemperature - Hypoxia      5.18     NA         NA     8   50.0  50.0
 9 HighTemperature - Hypoxia      5.18     NA         NA     9   55.3  55.3
10 HighTemperature - Hypoxia      5.18     NA         NA    10   52.6  52.6
# ℹ 21,590 more rows

norin.em |>
  group_by(contrast) |>
  median_hdci(Fit)

# A tibble: 3 × 7
  contrast                        Fit .lower .upper .width .point .interval
  <fct>                         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 HighTemperature - Hypoxia     50.9    28.9   73.2   0.95 median hdci     
2 HighTemperature - LowSalinity 41.8    22.2   62.6   0.95 median hdci     
3 Hypoxia - LowSalinity         -8.92  -30.0   11.6   0.95 median hdci

norin.em |>
  group_by(contrast, SMR_contr) |>
  median_hdci(Fit)

# A tibble: 9 × 8
  contrast                SMR_contr    Fit .lower .upper .width .point .interval
  <fct>                       <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 HighTemperature - Hypo…      3.93  64.3    53.7  76.3    0.95 median hdci     
2 HighTemperature - Hypo…      5.18  50.9    45.0  56.1    0.95 median hdci     
3 HighTemperature - Hypo…      6.43  37.0    25.8  49.3    0.95 median hdci     
4 HighTemperature - LowS…      3.93  49.4    33.1  66.9    0.95 median hdci     
5 HighTemperature - LowS…      5.18  41.7    33.3  50.4    0.95 median hdci     
6 HighTemperature - LowS…      6.43  34.3    16.5  51.7    0.95 median hdci     
7 Hypoxia - LowSalinity        3.93 -15.1   -33.6   4.15   0.95 median hdci     
8 Hypoxia - LowSalinity        5.18  -8.89  -18.6   1.33   0.95 median hdci     
9 Hypoxia - LowSalinity        6.43  -2.86  -23.2  16.0    0.95 median hdci

## norin.em |>
##     group_by(contrast) |>
##     summarize(P=sum(Fit>0)/n())
norin.em |>
  group_by(contrast, SMR_contr) |>
  summarise(
    P = mean(Fit > 0),
    P2 = 1 - P
  )

`summarise()` has grouped output by 'contrast'. You can override using the
`.groups` argument.

# A tibble: 9 × 4
# Groups:   contrast [3]
  contrast                      SMR_contr      P       P2
  <fct>                             <dbl>  <dbl>    <dbl>
1 HighTemperature - Hypoxia          3.93 1      0       
2 HighTemperature - Hypoxia          5.18 1      0       
3 HighTemperature - Hypoxia          6.43 1      0       
4 HighTemperature - LowSalinity      3.93 1      0       
5 HighTemperature - LowSalinity      5.18 1      0       
6 HighTemperature - LowSalinity      6.43 1.000  0.000417
7 Hypoxia - LowSalinity              3.93 0.0542 0.946   
8 Hypoxia - LowSalinity              5.18 0.0342 0.966   
9 Hypoxia - LowSalinity              6.43 0.397  0.603

10 Summary figures

brms

norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.brm4 |>
  emmeans(~ SMR_contr | TRIAL, at = norin.grid) |>
  as.data.frame()

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

head(newdata)

 SMR_contr TRIAL             emmean lower.HPD upper.HPD
  3.984319 HighTemperature 73.88891  63.57026  83.21900
  4.011756 HighTemperature 73.36914  63.20446  82.49482
  4.039192 HighTemperature 72.84601  62.88133  81.85916
  4.066629 HighTemperature 72.28818  63.18352  81.79213
  4.094065 HighTemperature 71.74851  62.75470  80.98760
  4.121502 HighTemperature 71.22935  62.35300  80.23893

Point estimate displayed: median 
HPD interval probability: 0.95

ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  theme_classic() +
  theme(
    legend.position = c(0.99, 0.99),
    legend.justification = c(1, 1)
  )

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

## The .fixed values are the predicted values without random effects
obs <- norin.brm4 |>
  augment() |>
  mutate(PartialObs = .fitted + .resid)

ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
  ## geom_point(data=norin,  aes(y=CHANGE), color='gray') +
  theme_classic()

g1 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
  theme_classic()
g2 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = norin, aes(y = CHANGE, color = TRIAL)) +
  theme_classic()
g1 + g2

11 References

norin <- norin |> mutate(
  FISHID = factor(FISHID),
  TRIAL = factor(TRIAL)
)

ggplot(norin, aes(y = CHANGE, x = TRIAL)) +
  geom_boxplot()
ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
  geom_smooth(method = "lm") +
  geom_point()
# ggplot(norin, aes(y=CHANGE, x=MASS, shape=TRIAL, color=TRIAL)) +
# geom_smooth(method='lm') + geom_point()
ggplot(norin, aes(y = CHANGE, x = as.numeric(FISHID), color = TRIAL)) +
  geom_point() +
  geom_line()

# ggplot(norin, aes(y=MASS, x=TRIAL)) + geom_boxplot()
ggplot(norin, aes(y = CHANGE, x = MASS, color = TRIAL)) +
  geom_point() +
  geom_smooth(method = "lm")

norin.rstanarm <- stan_glmer(CHANGE ~ (1 | FISHID) + TRIAL * SMR_contr + MASS,
  data = norin,
  prior_PD = TRUE,
  iter = 5000, warmup = 2000, chains = 3, thin = 5, refresh = 0
)
prior_summary(norin.rstanarm)

posterior_vs_prior(norin.rstanarm,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y"), pars = c("(Intercept)")
)
ggpredict(norin.rstanarm, ~ TRIAL * SMR_contr) |> plot(show_data = TRUE)

norin.rstanarm |> get_variables()
plot(norin.rstanarm, "mcmc_trace", regex_pars = "^.Intercept|TRIAL|SMR|MASS|[sS]igma")
plot(norin.rstanarm, "mcmc_acf_bar", regex_pars = "^.Intercept|TRIAL|SMR|MASS|[sS]igma")
plot(norin.rstanarm, "mcmc_rhat_hist", regex_pars = "^.Intercept|TRIAL|SMR|MASS|[sS]igma")
plot(norin.rstanarm, "mcmc_neff_hist", regex_pars = "^.Intercept|TRIAL|SMR|MASS|[sS]igma")

# norin.rstan1 = stan_glmer(CHANGE ~ (TRIAL|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
#                          iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.rstanarm1 <- stan_glmer(CHANGE ~ (SMR_contr | FISHID) + TRIAL * SMR_contr + MASS,
  data = norin,
  prior_PD = FALSE,
  iter = 5000, warmup = 2000, chains = 3, thin = 5, refresh = 0, cores = 3
)
norin.rstanarm1 <- update(norin.rstanarm1, prior_PD = FALSE)



norin.rstanarm2 <- stan_glmer(CHANGE ~ (TRIAL * SMR_contr | FISHID) + TRIAL * SMR_contr + MASS,
  data = norin,
  prior_PD = FALSE,
  iter = 5000, warmup = 2000, chains = 3, thin = 5, refresh = 0, cores = 3
)

posterior_vs_prior(norin.rstanarm1,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y"), pars = c("(Intercept)")
)

ggpredict(norin.rstanarm1, ~ TRIAL * SMR_contr) |> plot(show_data = TRUE)

norin.rstanarm1 |> get_variables()
plot(norin.rstanarm1, "mcmc_trace", regex_pars = "^.Intercept|TRIAL|^SMR|MASS|[sS]igma")
plot(norin.rstanarm1, "mcmc_acf_bar", regex_pars = "^.Intercept|TRIAL|^SMR|MASS|[sS]igma")
plot(norin.rstanarm1, "mcmc_rhat_hist", regex_pars = "^.Intercept|TRIAL|^SMR|MASS|[sS]igma")
plot(norin.rstanarm1, "mcmc_neff_hist", regex_pars = "^.Intercept|TRIAL|^SMR|MASS|[sS]igma")


(l.1 <- loo(norin.rstanarm))
(l.2 <- loo(norin.rstanarm1))
loo_compare(l.1, l.2)


preds <- posterior_predict(norin.rstanarm, nsamples = 250, summary = FALSE)
norin.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = norin$CHANGE,
  fittedPredictedResponse = apply(preds, 2, median)
)
plot(norin.resids)


g <- ggpredict(norin.rstanarm) |> plot()
do.call("grid.arrange", g)

# ggemmeans(norin.rstan, ~TRIAL)

summary(norin.rstanarm)
nms <- norin.rstanarm1 |> get_variables()
wch <- grep("^.Intercept|TRIAL|^SMR|[sS]igma", nms)
tidyMCMC(norin.rstanarm$stanfit,
  conf.int = TRUE, conf.method = "HPDinterval",
  rhat = TRUE, ess = TRUE, pars = nms[wch], estimate.method = "median"
)

tidyMCMC(norin.rstanarm1$stanfit,
  conf.int = TRUE, conf.method = "HPDinterval",
  rhat = TRUE, ess = TRUE, pars = nms[wch], estimate.method = "median"
)


norin.grid <- with(norin, list(SMR_contr = seq(min(SMR_contr), max(SMR_contr), len = 100)))
newdata <- emmeans(norin.rstanarm, ~ TRIAL | SMR_contr, at = norin.grid) |> as.data.frame()
head(newdata)
ggplot(newdata, aes(y = emmean, x = SMR_contr, color = TRIAL)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3, color = NA) +
  geom_line()

norin.grid <- with(norin, list(SMR_contr = c(min(SMR_contr), mean(SMR_contr), max(SMR_contr))))

emmeans(norin.rstan, pairwise ~ TRIAL | SMR_contr, at = norin.grid)

norin.em <- emmeans(norin.rstan, pairwise ~ TRIAL | SMR_contr, at = norin.grid)$contrast |>
  gather_emmeans_draws() |>
  mutate(Fit = .value)
norin.em
norin.em |>
  group_by(contrast) |>
  median_hdci(Fit)
norin.em |>
  group_by(contrast, SMR_contr) |>
  median_hdci(Fit)
## norin.em |>
##     group_by(contrast) |>
##     summarize(P=sum(Fit>0)/n())
norin.em |>
  group_by(contrast, SMR_contr) |>
  summarize(P = mean(Fit > 0))


bayes_R2(norin.rstanarm, re.form = NA) |> median_hdi()
bayes_R2(norin.rstanarm, re.form = ~ (1 | FISHID)) |> median_hdi()
# bayes_R2(norin.rstan1, re.form=~(SMR_contr|FISHID)) |> median_hdi

11.1 brms

norin <- norin |> mutate(
  FISHID = factor(FISHID),
  TRIAL = factor(TRIAL)
)

ggplot(norin, aes(y = CHANGE, x = TRIAL)) +
  geom_boxplot()
ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
  geom_smooth(method = "lm") +
  geom_point()
ggplot(norin, aes(y = CHANGE, x = MASS, shape = TRIAL, color = TRIAL)) +
  geom_smooth(method = "lm") +
  geom_point()
ggplot(norin, aes(y = CHANGE, x = as.numeric(FISHID), color = TRIAL)) +
  geom_point() +
  geom_line()

## ggplot(norin, aes(y=MASS, x=TRIAL)) + geom_boxplot()
## ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) + geom_point() + geom_smooth(method='lm')

norin |>
  group_by(TRIAL) |>
  summarise(
    median(CHANGE),
    mad(CHANGE)
  )
priors <- prior(normal(50, 20), class = "Intercept") +
  prior(normal(0, 10), class = "b") +
  prior(gamma(2, 1), class = "sigma") +
  prior(gamma(2, 1), class = "sd")

norin.form <- bf(CHANGE ~ (1 | FISHID) + TRIAL * SMR_contr + MASS,
  family = gaussian
)

norin.brm1 <- brm(norin.form,
  data = norin,
  prior = priors,
  sample_prior = "yes",
  iter = 5000, warmup = 2000,
  chains = 3, thin = 5, refresh = 0
)

norin.brm1 |> get_variables()
pars <- norin.brm1 |> get_variables()
wch <- grepl("^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd", pars, perl = TRUE)

stan_trace(norin.brm1$fit, pars = pars[wch])
stan_ac(norin.brm1$fit, pars = pars[wch])
stan_rhat(norin.brm1$fit, pars = pars[wch])
stan_ess(norin.brm1$fit, pars = pars[wch])

## mcmc_plot(norin.brms,  type='trace',
##          regex_pars='^b.*|sigma|^sd')
## mcmc_plot(norin.brms,  type='trace',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brms,  type='acf_bar',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brms,  type='rhat_hist',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brms,  type='neff_hist',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')

preds <- posterior_predict(norin.brm1, nsamples = 250, summary = FALSE)
norin.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = norin$CHANGE,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = FALSE
)
plot(norin.resids)
# norin.rstan1 = stan_glmer(CHANGE ~ (TRIAL|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
#                          iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.form <- bf(CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + MASS,
  family = gaussian
)
norin.brm2 <- brm(norin.form,
  data = norin,
  prior = priors,
  sample_prior = "yes",
  iter = 5000, warmup = 2000,
  chains = 3, thin = 5, refresh = 0, cores = 3,
  control = list(adapt_delta = 0.99)
)

norin.brm2 |> get_variables()

pars <- norin.brm2 |> get_variables()
## wch <- grepl('^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd', pars, perl=TRUE)
wch <- grepl("^b_.*|[sS]igma|^sd_.*", pars, perl = TRUE)

stan_trace(norin.brm2$fit, pars = pars[wch])
stan_ac(norin.brm2$fit, pars = pars[wch])
stan_rhat(norin.brm2$fit) # , pars=pars[wch])
stan_ess(norin.brm2$fit) # , pars=pars[wch])
## mcmc_plot(norin.brm2,  type='trace',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brm2,  type='trace',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brm2,  type='acf_bar',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brm2,  type='rhat_hist',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brm2,  type='neff_hist',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')

(l.1 <- loo(norin.brm1))
(l.2 <- loo(norin.brm2))
loo_compare(l.1, l.2)


preds <- posterior_predict(norin.brm2, nsamples = 250, summary = FALSE)
norin.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = norin$CHANGE,
  fittedPredictedResponse = apply(preds, 2, median)
)
plot(norin.resids)

g <- norin.brm2 |>
  conditional_effects() |>
  plot(points = TRUE, ask = FALSE)
library(patchwork)
g[[1]] + g[[2]] + g[[3]] + g[[4]]


## g=ggpredict(norin.brms1) |> plot
## library(patchwork)
## g[[1]] + g[[2]] + g[[3]]

## do.call('grid.arrange', g)

ggemmeans(norin.brm2, ~TRIAL) |> plot()

summary(norin.brm2)

tidyMCMC(norin.brm2$fit,
  conf.int = TRUE, conf.method = "HPDinterval",
  rhat = TRUE, ess = TRUE, estimate.method = "median"
) |>
  slice(1:11)

pars <- norin.brm2 |> get_variables()
wch <- grep("^b.Intercept|TRIAL|^b.*SMR|[sS]igma|^sd", pars)
tidyMCMC(norin.brms1$fit,
  conf.int = TRUE, conf.method = "HPDinterval",
  rhat = TRUE, ess = TRUE, pars = pars[wch], estimate.method = "median"
)

bayes_R2(norin.brm2, re.form = NA, summary = FALSE) |>
  median_hdci()
bayes_R2(norin.brm2, re.form = ~ (1 | FISHID), summary = FALSE) |>
  median_hdci()
bayes_R2(norin.brm2, re.form = ~ (TRIAL | FISHID), summary = FALSE) |>
  median_hdci()

emmeans(norin.brm2, pairwise ~ TRIAL)


norin.em <- norin.brm2 |>
  emmeans(~TRIAL) |>
  pairs() |>
  gather_emmeans_draws() |>
  mutate(Fit = .value)

norin.em |>
  group_by(contrast) |>
  median_hdi()

norin.em |>
  ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  stat_slab(aes(
    x = .value, y = contrast,
    fill = stat(ggdist::cut_cdf_qi(cdf,
      .width = c(0.5, 0.8, 0.95),
      labels = scales::percent_format()
    ))
  ), color = "black") +
  scale_fill_brewer("Interval", direction = -1, na.translate = FALSE) +
  theme_bw()

norin.em |>
  group_by(contrast) |>
  summarize(P = mean(Fit > 0))


norin.grid <- with(norin, list(SMR_contr = c(
  min(SMR_contr),
  mean(SMR_contr),
  max(SMR_contr)
)))

norin.em <- norin.brm2 |>
  emmeans(~ TRIAL | SMR_contr, at = norin.grid) |>
  pairs() |>
  gather_emmeans_draws()

norin.em |> head()
norin.em |>
  group_by(contrast, SMR_contr) |>
  median_hdi()

norin.em |>
  group_by(contrast, SMR_contr) |>
  summarize(P = mean(.value > 0))

norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.brm2 |>
  emmeans(~ SMR_contr | TRIAL, at = norin.grid) |>
  as.data.frame()
head(newdata)
partial.obs <- norin |>
  mutate(
    Pred = predict(norin.brm2, re.form = NA, summary = TRUE)[, "Estimate"],
    Resid = resid(norin.brm2)[, "Estimate"],
    Obs = Pred + Resid
  )
ggplot(newdata, aes(y = emmean, x = SMR_contr, color = TRIAL)) +
  geom_point(data = partial.obs, aes(y = Obs)) +
  ## geom_point(data=partial.obs, aes(y=CHANGE), shape=2) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3, color = NA) +
  geom_line()

References

Norin, Tommy, Hans Malte, and Timothy D. Clark. 2015. “Differential Plasticity of Metabolic Rate Phenotypes in a Tropical Fish Facing Environmental Change.” Functional Ecology 30 (3): 369–78. https://doi.org/10.1111/1365-2435.12503.